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 frequency 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 analyzing 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 relation of the PRC shape, voltage time course, and stability of the resultant synchronous and antisynchronous solutions.
The theory of weakly coupled oscillators (Winfree, 1967; Neu, 1979a, 1979b; Kuramoto, 1984; Ermentrout & Kopell, 1986, 1991; Hoppensteadt & Izhikevich, 1997; Brown, Moehlis, & Holmes, 2004) provides a mechanism to relate the phase response curve (PRC) of a neuron and its voltage to the stability of the emergent phase-locked states such as synchrony of a network. The PRC and the voltage are not totally independent of one another (Ermentrout & Kopell, 1986, 1991; Hoppensteadt & Izhikevich, 1997; Brown et al., 2004), but at the same time, their shapes are not easily predictable from each other. These shapes could be affected by a steady stimulus, drug application, or modulation of ion channel conductances (Ermentrout, Pascal, & Gutkin, 2001; Netoff et al., 2005; Tateno & Robinson, 2007). For example, causing the frequency of a neuron to change would also have caused a number of other changes in the shapes of the PRC and the voltage: a change of the spike width, spike height, spike maximum hyperpolarization, spike threshold, a possible change of the slope of the depolarizing phase, change of the PRC skewness, and a change of the maximum delay and advancement of the PRC. A detailed and extensive understanding of which of these parameters is more sensitive in changing the dynamical behavior of a network is generally lacking. Considerable progress is made when the neurons are synaptically coupled and are posited near the bifurcation point (Van Vreeswijk, Abbott, & Ermentrout, 1994; Hansel, Mato, & Meunier, 1995; Ermentrout, 1996) where the voltage time course does not significantly affect the stability of synchrony and antisynchrony, leading to a somewhat general theory that relates the PRC shapes to the network behavior.
However, when the neurons are coupled by gap junctions or electrical connections, the task becomes more difficult. In fact, very few general rules that relate the PRC and voltage shapes to the network behavior are available in the literature. Most studies have focused on leaky integrate-and-fire and related models (Chow & Kopell, 2000; Lewis & Rinzel, 2003), and the predictions are confined to specific forms of the PRCs corresponding to those models. This is still an unsolved problem requiring a detailed study. Lewis and Skinner (2012) summarize the importance of this problem by stating that in order to obtain detailed insights into the dynamical and biophysical mechanisms underlying network behavior, we must “determine how the shapes of Z and VLC affect phase-locking, i.e., study how the shapes of the functions Z and VLC combine to influence the phase-locking states.” (Z and VLC are, respectively, the PRC and the voltage profile of the neuron between two consecutive spikes.)
Cell-to-cell gap-junction-mediated electrical coupling is important because it is not only found in a number of brain neurons, including striatal fast-spiking interneurons (Koos & Tepper, 1999; Klaus et al., 2011), neocortex (Galarreta & Hestrin, 1999; Gibson, Beierlein, & Connors, 1999; Mancilla, Lewis, Pinto, Rinzel, & Connors, 2007), thalamic reticular cells (Landisman et al., 2002), thalamic relay neurons during development (Lee, Cruikshank, & Connors, 2010), and CA1 hippocampal neurons that could generate synchronous discharges characteristic of seizures (Valiante, Perez Velazquez, Jahromi, & Carlen, 1995; Carlen et al., 2000) but is also common in several other biological preparations such as sinoatrial node cells (Jalife, 1984; Verheijck et al., 1998; Demir, Clark, & Giles, 1999). When the neurons are in a oscillatory state either autonomously or in response to external stimuli, a synchronous state may emerge in the connected network (Kepler, Marder, & Abbott, 1990; Traub et al., 2003; Fuentealba et al., 2004; Hestrin & Galarreta, 2005; Mancilla et al., 2007 (also see simulations in Gao & Holmes, 2007; Ostojic, Brunel, & Hakim, 2009). Electrical coupling has also been found to result in a failure of synchrony (Bou-Flores & Berger, 2001). A number of models have been proposed to explain the mechanism of synchrony or its failure (Chow & Kopell, 2000; Lewis & Rinzel, 2003; Nomura, Fukai, & Aoyagi, 2003; Bem, Le Feuvre, Rinzel, & Meyrand, 2005; Pfeuty, Mato, Golomb, & Hansel, 2003). When neurons are coupled electrically, the coupling lasts the duration of the voltage time course. Although this dependence of the coupling on the voltage shape makes the task of obtaining general rules more difficult, some earlier studies succeeded in the parameterization of spike width, spike height, and frequency while keeping the rest of the spike profile and the PRC shape intact (Chow & Kopell, 2000; Lewis & Rinzel, 2003).
In this study, we are interested in the stability of phase-locked states of two coupled neurons that oscillate at identical natural frequency in their uncoupled state. The phase-locked states we are interested in are synchrony and antisynchrony. The synchronous state is characterized by the oscillating neurons maintaining no spike time difference, whereas the antisynchrony is characterized by the neurons maintaining a spike time difference equal to half their oscillation period. The model of the neuron, however, is not based on a dynamical evolution of independent variables. It is represented by its PRC and its voltage shape, which together with the knowledge of the coupling mechanism are sufficient to predict the emerging phase-locked states and their linear stability. Because we aim to obtain some general rules relating the PRC and the voltage shapes to the phase-locked states, we parameterize the shapes of both the PRC and the voltage, which is facilitated by not imposing any functional relationship between them in the analysis. This dissociation of the shapes would extend the applicability of the results to a wider class of PRCs and voltage shapes that are found in both experiments and models. This also gives the ability to independently evaluate the sensitivity of the phase-locked states to changes in the parameters that define the shapes of the PRC and the voltage. But we first need to make a reasonable choice of these shapes.
Phase response curves are usually classified as type 1 or type 2 (Hansel et al., 1995) depending on whether they show only phase advancement or that in addition to phase delay. The classic Hodgkin-Huxley neuron model, when driven using a steady current to oscillate, displays a type 2 PRC. We aim to capture not only type 2 PRCs of this kind but also other type 1 PRCs that may be obtained by altering the shape parameters. Fitting experimental PRCs with Fourier modes (Galán, Ermentrout, & Urban, 2005; Tsubo, Teramae, & Fukai, 2007) or polynomials (Netoff et al., 2005) is a common practice in neuroscience, so it would seem an obvious choice to formulate shapes comprising Fourier modes or polynomials. However, the functional form of these shapes is less crucial than the ability to parameterize those shapes to span a wide range of PRCs. Considering that at least three Fourier modes are recommended for fitting most experimental PRCs, the number of parameters required to achieve a wide range of both type 1 and type 2 PRCs can be large. Similarly a third or higher-order polynomial could be considered, but it still suffers from the same disadvantage as that of the Fourier series.
Instead we use a piecewise linear functional form to model the PRCs, requiring effectively only two independent parameters: a skewness parameter that determines how far the maximum phase advancement is from the spike initial phase and a type parameter that determines the relative magnitude of the maximum phase delay with respect to the maximum phase advancement. The skewness parameter could range from zero, when the maximum phase advancement of the PRC is at half-oscillation period, to the oscillation period itself, when the PRC has its peak near the oscillation period. The type parameter could be zero or positive, making it a type 1, or negative, making it a type 2. As will be demonstrated, this choice achieves our goals of spanning a wide range of shapes while being still accurate in predicting the nature of stability of the phase-locked states. Voltage profiles broadly comprise spike downstroke, the depolarization phase, and the spike upstroke. We again consider piecewise linear functional forms to model the voltage time course that is shaped by few parameters: the normalized spike width parameter that could range from 0 to its maximum allowed value and the three amplitude parameters—spike width, spike height, and spike threshold. All of these could be altered freely. Note that though linear in individual segments, the complete PRC and voltage shapes themselves are nonlinear. The PRC and the voltage shape parameters are fully explored in the permissible ranges of parameters to delineate the regions of stability of both synchrony and antisynchrony.
In section 2, the PRC and the voltage shapes are formulated, and the method to determine the stability of the phase-locked states is described. The stability of synchronous and antisynchronous states is studied in section 3, when the spike width is zero, and in section 4, when the spike width is nonzero. The results are discussed in relation to previous studies and are analyzed in section 5. We obtain analytical results and provide explicit relations for the boundaries of the phase-locked states. These boundaries are rigorously verified when feasible by examining the corresponding eigenvalue transitions. The results are summarized in Figure 5 for zero spike width and in Figures 10 and 11 for nonzero spike width.
2. Model and Methods
2.1. Model for the Voltage Time Course
2.2. Model for the Phase Response Curve
A popular method of measuring phase response curves experimentally is by computing the relative shifts of the spike times in response to brief current inputs placed during the time course in between the spike times. However, such current inputs also cause similar phase shifts in other independent variables that define the spiking dynamics of the neuron. The profiles of these phase shifts are different in different variables and together form different components of the adjoint solution that can be computed numerically for any spiking neuron model (Izhikevich, 2007). This is also referred to as the phase response curve, linear response function, or sensitivity function (Winfree, 1967). However, since the electrical coupling depends on only the voltage, the voltage component of the phase response curve [Z(t)] (i.e., the phase shifts measured in the voltage variable) is of interest to us, and we formulate it by five piecewise linear curves (see Figure 1). It is parameterized by three parameters A, B, and C. Of these three, A is a time parameter, and B and C (which can indeed be combined into a single parameter B/C without loss of generality) define the magnitude levels of the PRC. In addition to these three parameters, the time period T controls both the voltage and the PRC. Spike rise time W/2 is also used in specifying the PRC.
PRC with large A is said to have large skewness, and that with no or zero skewness is characterized by A=0. Many experimental PRCs show the maximum phase advancement tilted more toward higher values of the phases, that is, they have large skewness. Mimicking the HH model that has a very small response near early phases, we have assumed the PRC to be zero from spike peak time to a time half of the maximum delay time of the PRC. Also, the PRC is assumed to be zero during the rise time of the spike. The zero regime during the action potential is essentially a depiction of the fact that the identical phase lines (i.e., isochrons) drawn in the phase plane of voltage versus any other variable are parallel to the voltage axis; thus, perturbations imparted to the voltage (which are used in determining the PRC) cause negligible phase shifts.
2.3. Stability of Synchrony and Antisynchrony
The regions and the relative placement of V and Z in some of these regions are illustrated in Figure 2B. We have detailed the procedure and the eigenvalue expressions for all the components in, respectively, appendix B and appendix E. An examination of the relative placement of V and Z can give us an estimate of the sensitivity of the synchronous and antisynchronous states to the parameters that define the PRC and voltage shapes.
In the following two sections, we analyze these eigenvalues to determine where synchronous and antisynchronous states are located in the parameter regimes. Explicit curves for the boundaries of these states are derived. When possible, the derivation is substantiated by examining the actual transitions of the eigenvalues.
3. Zero Spike Width
In this section, conditions for synchrony and antisynchrony are derived analytically for the case of zero spike width (W=0). An example set of type 1 PRCs is displayed in Figure 3b for a range of PRC skewness levels. These are obtained by setting B=0. The parameter region is represented by region c in the (W,T) plane of Figure 2A and is illustrated by an example in case c, type I. We clearly see from the lack of any shaded region there (see appendix A for a description) that no eigenvalue component is positive for any level of skewness, and the nonzero eigenvalue components and , displayed in Figure 3c, and hence the total eigenvalue are negative and thus the synchrony is stable for all A. For the antisynchronous state we refer to regions b and e of the (W,T) plane of Figure 2B, and the corresponding illustrations in cases b and e of type I. There we see that the only eigenvalue component that is positive, and hence contributes to instability of antisynchrony, is due to the spike downstroke, ( becomes zero at B=0). We will see that this is the biggest in magnitude among all the components due to the contribution of the downstroke with infinite slope in the limit of zero spike width. All of the other components are negative. These are illustrated in Figure 3d. As the skewness is increased, the discontinuity in the voltage occurs at smaller and smaller levels of PRC, leading to a diminishing effect of it, and thus the total eigenvalue becomes negative beyond certain level of skewness, causing the antisynchrony to become stable. A numerical bifurcation diagram obtained by solving for the slopes of directly from equation 2.4 (see Figure 3e) confirms these observations. The stable synchrony and antisynchrony regions as a function of the skewness are shown in Figure 3f and will be rigorously shown in this section.
An example of type 2 PRCs for a range of A is shown in Figure 4b. B/C is set to −0.5 for illustration. For the eigenvalue components contributing to the stability of synchrony, we refer to case c, type II, in Figure 2A. There we see that the negative lobe of the PRC can potentially destabilize synchrony (see appendix B for a description) if its contribution surpasses in magnitude that of the other PRC regions. Of the two shaded segments, the Z2(t) segment of the PRC results in a positive eigenvalue component () for any level of skewness, and its magnitude increases linearly with skewness. At the same time, the other PRC segments span a smaller range along the phase, and thus their contribution to the negative eigenvalue diminishes (see Figure 4c), leading to the total eigenvalue becoming positive above a critical skewness and causing the synchrony to become unstable. For the stability of antisynchrony, we refer to cases b and e of type II in Figure 2B. At small skewness, the spike discontinuity occurs at a positive level of PRC, and thus it contributes a positive component to the eigenvalue (sum of due to upstroke and due to downstroke) that overpowers the other stabilizing components of the eigenvalue. This effectively makes the antisynchrony unstable. But as the skewness is increased such that the discontinuity occurs closer to the zero crossing of the PRC, its effect in destabilizing the antisynchrony diminishes, leading to a stable antisynchrony. At extremely large skewness, the discontinuity occurs at a negative PRC level when the stabilizing roles of the spike upstroke and downstroke are reversed (represented now by and ). But at the same time, the rising phase of the voltage convolved with the negative PRC regime ( and ) adds to the positive level of the eigenvalue. Consequently at very large skewness, the antisynchrony could become unstable again (see Figure 4d). A numerical bifurcation diagram as a function of skewness (see Figure 4e) confirms that the synchrony becomes unstable at large skewness, whereas antisynchrony is unstable for small A/T but becomes stable for large A/T before becoming unstable again at extremely large A/T. The analytical boundaries for synchrony () and the antisynchrony ( and ) that are derived later in this section are depicted in Figure 4f. We in fact derive the stability regions in the plane of skewness and type parameter; These are depicted in Figure 5a. The type 1 PRCs always show stable synchrony, and type 2 PRCs are capable of displaying stable synchrony at all skewness levels. Both type 1 and type 2 PRCs exhibit stable antisynchrony at large skewness, with type 2 PRCs showing instability at very large skewness levels. Type 2 PRCs may also lose stability at very large skewness. These are verified by numerically computed in Figures 5b to 5d. We next present the analytical stability boundaries.
Two stability criteria will emerge for antisynchrony, corresponding to the case of small skewness (see case b, Figure 2B) and the case of large skewness (see case e, Figure 2B). In the first case, , and the eigenvalue determines the stability, and in the second case, , where determines the stability. We will find that there are four segments that form the boundary of stable antisynchrony. Type 1 PRCs at moderately large skewness () may lead to unstable antisynchrony. They can also lead to stable antisynchrony for very small skewness () provided the type parameter is large positive. Furthermore, we will see that at very large skewness (), a type 2 PRC can lead to unstable antisynchrony. There are two disjoint regions of stability in (A, B) parameter space.
4. Nonzero Spike Width
In this section, we investigate the effect of spike width and the dependence of stability criteria on it. The effects of the spike width in type 1 PRCs is illustrated in Figure 6 (compare with Figure 3). The parameter W/T is set at 0.15 in the figure. For studying the stability of synchrony, in the limit of zero spike width, we could consider the parameter point as belonging to region c as we did earlier or to region a. But for finite spike width, we must consider regions a and d in the (W, A) space of Figure 2A and the corresponding illustrations in cases a and d. When A/T=0, the PRC is nonzero near early phases, and thus the spike downstroke could cause instability of synchrony (see Figure 6c). As the skewness is increased, the spike downstroke effect diminishes because of the either small or zero level of the PRC, and consequently the synchrony acquires stability. Again in contrast to the case of zero spike width, the antisynchrony can acquire stability even at A/T=0. At finite spike width, the spike downstroke, case a, type I in Figure 2B, that earlier destabilized antisynchrony could become less effective ( in Figure 6d) if it extends to the region of the PRC that has negative slope. Consequently the total eigenvalue can become slightly negative for certain W/T, leading to a stable antisynchrony. But as the skewness is increased, the spike downstroke that contributes to a positive eigenvalue occurs at PRC levels that are larger than that during the spike upstroke, leading to unstable antisynchrony. At extremely large spike width, stability can return as quantified by the eigenvalue components in Figure 6d. A numerical one-parameter bifurcation diagram as a function of A/T obtained by computing the from equation 2.4 (see Figure 6e) confirms these observations. The corresponding boundaries of stability are shown in Figure 6f. In the illustration we have set B/C=0. For finite positive B/C, as will be seen from the expressions derived later, the movement of the boundaries are as indicated in the side panel of Figure 6f.
The effect of spike width on type 2 PRCs is illustrated in Figure 7. W/T is again set to 0.15 and B/C=−0.5 (see Figures 7a and 7b). At these parameter values, the stability of synchrony is similar to that at zero spike width. The spike downstroke could still cause instability of synchrony ( in case a, type II of Figure 2A) even at small skewness, but due to the narrow region between the zero crossing of the PRC and the time of the maximum hyperpolarization that contributes to the positive part of , it would require a large negative B/C to cause instability. But in regions b and c, part or all of the negative lobe of the PRC (see Figure 2A, cases b and c, type II) contributes to a positive eigenvalue (see Figure 7c) causing the synchrony to become unstable. The effect of spike width on the antisynchrony could be more drastic than on synchrony. The antisynchrony becomes unstable for most of the skewness except at large and very small skewness. The components and (see Figure 7d) (and their corresponding segments in regions f, d, g, and e) together dominate the contribution to the total eigenvalue in causing the instability. A numerical one-parameter bifurcation diagram in Figure 7e as a function of skewness confirms these observations. The stability boundaries and their dependence on the level of B/C are depicted in Figure 7f.
Before we present the analytical results, we illustrate the role of spike threshold together with spike width at a small (see Figure 8, A/T=0.2) and a large (see Figure 9, A/T=0.6) value of skewness. Comparing Figures 3e and 6e, we see that finite spike width could cause an instability of synchrony at small skewness. From these figures we also notice that the range of stable antisynchrony appearing near large skewness levels is moved further to larger skewness. The movement of eigenvalues for antisynchrony is shown in Figure 8a, and the corresponding bifurcation diagram in the space of W/T and a3/a2 is shown in Figure 8b. When the synchrony is unstable, increasing the spike threshold such that a3/a2 crosses a critical level causes the synchrony to become stable because the destabilizing effect of downstroke is countered by the contribution of the stabilizing segments of the PRC that are now amplified by a larger slope of the voltage segment. computed at a parameter point where both the synchrony and antisynchrony are unstable is shown in Figure 8c, and a one-parameter bifurcation diagram, Figure 8d, verifies the stability diagram. A stability diagram at the same level of skewness for a type 2 PRC (here, by setting B/C=−0.5) is shown in Figure 8e and is numerically verified in Figure 8f.
Large skewness causes the spike downstroke to have less effect on synchrony and thus is generally favorable to stable synchrony. Except when the spike threshold is small, such that the slope of the depolarizing phase is very small, we find a stable synchrony (see Figure 9a) for all values of W/T. Large skewness coupled with small or zero spike width also makes the spike profile itself be less effective (see cases d and e in Figure 2B) in the antisynchrony. The spike downstroke contributes to instability in type 1 PRCs but is now diminished and contributes to stability in type 2 PRCs. Thus, we may expect a stable antisynchrony at a small spike width. It becomes unstable at large spike width and frequency (see Figure 9a). A numerical bifurcation diagram in Figure 9b verifies this. The eigenvalue components illustrated for a type 1 PRC (see Figure 9c) for the antisynchrony clearly reveal that the stabilizing segments (see and ) at small spike width or frequency are not countered by the spike contribution. A numerical computation of (see Figure 9d) reveals that the antisynchrony is stable at small W/T and becomes unstable at large W/T, while the synchrony is stable at all W/T.
The four cases (a–d) depicted in the (W, A) space in Figure 2A are treated here, and the corresponding regions of stability for synchrony are derived. Since each case constitutes only a portion of the parameter space, the critical curves in the space of (W, A) or (A, B) are composed of one or more of these critical curves. The eigenvalues are linear in the type parameter, B, but even then, the expressions are not as simple as those in the previous section, and the eigenvalue transition conditions cannot be easily ascertained except in the case of very large skewness (see case c). In the other three cases, we provide the critical curves by solving the critical eigenvalue equations and list out expressions for the eigenvalue transitions. The regions may be verified by directly plotting the eigenvalue and checking for regimes where it is negative. The results from this section are presented in Figures 6 and 7, which also display the explicit movement of eigenvalues as one of the parameters is varied, and in Figures 8 to 11 that are plotted after verifying the eigenvalues (not displayed).
To illustrate the effect of finite spike width on the phase-locked states, for type 1 (see Figure 6b) and type 2 (see Figure 7b) PRCs, we compute numerically a one-parameter bifurcation diagram as a function of PRC skewness for B/C=0 (see Figure 6e) and for B/C=−0.5 (see Figure 7e) at W/T=0.15. Type 1 PRCs showed a loss of synchrony at small skewness when the effect of spike downstroke dominates (), and type 2 PRCs showed a loss of synchrony at large skewness where the effect of the negative PRC lobe dominates (). The stability of synchrony in these diagrams is also verified by plotting the eigenvalues, respectively, in Figures 6c and 7c. Corresponding stability regions as a function of PRC skewness are shown in Figures 6f and 7f. The boundaries for type 1 neurons are defined by (see case a below) and (see case b), and for type 2 neurons, the stability boundaries are defined by , and (see case d). The stability boundaries for these four cases are given below. In case c where , the stability regions are derived exactly, but in the other three cases, the regions are not derived; only the boundary curves are listed due to the inherent complexity in the analysis, but we show the stability regions in (A/T, B/C) and (W/T, B/T) parameter spaces in Figures 10 and 11 based on the eigenvalue movements.
4.1.1. Very Large Skewness (A≥ 4W), Case c
4.1.2. Large Skewness (2W≤ A> 4W), Case b
4.1.3. Intermediate Skewness (4W-T ≤ A < 2W), Case a
4.1.4. Small Skewness (0 ≤ A < 4W-T), Case d
The curves , , and are illustrated in Figure 10 for select values of W/T at a3/a2=0.2234, and in Figure 11 for two skewness levels for three values of a3/a2. They are also illustrated in Figures 8 and 9. Figures 8 and 9 essentially indicate that the stable synchronous region expands in size to larger W/T with increasing skewness. This fact is also corroborated by the two-parameter stability regions in Figure 11. Increasing the type parameter B/C to negative values also increases this range, but beyond some level (, which may be beyond realistic B/C levels in models or experiments), synchrony loses stability and is largely confined to higher values of W/T.
It is difficult to derive conditions for the stability of antisynchrony in the general case due to the inherent complexity of the coefficients that arise in the stability conditions, but the critical curves can be written easily by solving corresponding eigenvalue equations. We list them here without rigorous proof of the stable regions. The type parameter (B) appears linearly in all the eigenvalue equations, and thus the critical curves expressed in terms of B are simpler without multiplicity than those expressed in terms of A or W. The critical curves are expressed in each of the 17 cases (case a, , case q), which are all marked in the (W, A) space in Figure 2. The arrangement of Z(t) and V(t) for some of the cases is also displayed. The main results are summarized in the two-parameter stability regions in Figures 10 and 11.
For small skewness (A<(T−W)/2 in (W, A) plot in Figure 2B), six cases could form boundaries of antisynchrony: cases b, a, j, k, l, and m. For example, four of these cases contribute to the boundary at A/T=0.2 (see Figures 8b and 8e). The combined value of the components ( and ) due to the spike downstroke for A>0 is comparable to that () for A=0 but could be higher because part of the destabilizing downstroke occurs in the ramp-up phase of the PRC, and thus the stability of antisynchrony is delayed until larger W/T (see Figure 8b). A similar effect is seen for negative B (see Figure 8e), as well as in Figure 11 where the stability boundary is defined by , , , , , and . In each of these cases, the eigenvalue components are given in appendix E, and using the eigenvalue equation 2.17, the critical curves are obtained by solving the equations , , . The critical curves written in terms of normalized type parameter are given in equations C.1 when .
The critical curve at A/T=0.2, B/C=0.2 (see Figure 8b) can be approximated by the first term in the Taylor series: 4.9(W/T−0.165). Inverting this term, we obtain an approximate stability criterion as W/T>0.165+0.2a3/a2. Similarly, the critical curve at A/T=0.2, B/C=−0.25 (see Figure 8e) can be approximated by the first term in the Taylor series: 4(W−0.22). Again inverting this term, we get the approximate stability criterion as W/T>0.22+0.25a3/a2. The stability of antisynchrony at high frequencies is verified numerically in Figures 8d and 8f.
Finite spike width and frequency helped in stabilizing antisynchronous state even in the absence of skewness (see Figure 6) essentially due to the sharp spike upstroke (). But the level of the PRC during the upstroke of V(t−T/2) decreases if the PRC acquires a finite skewness, causing the negative eigenvalue component to reduce its magnitude (see Figure 7d) resulting in a loss of antisynchronous state (see Figures 7e and 7f). But this loss can be delayed until larger values of A/T if we prevent the reduction of the PRC due to skewness, which can be achieved by increasing the B/C value as depicted in the side panel of Figure 7f. Reducing the B/C to negative values has the opposite effect (see Figure 7d and 7f) because the PRC level during the upstroke is further reduced.
For intermediate levels of skewness, the stability is defined by eigenvalues in the cases c, f, i, n, and o (W/T, A/T) space in Figure 2B. Similar to the curves derived above, we solve the equations , to obtain curves on which the eigenvalues are zero. Thus, the critical curves in terms of normalized type parameter B in the case when are given in equations C.2. And finally for large skewness (A>T/2), the critical curves are derived from cases e, d, g, h, p, and q. Following the same procedure as described above, we arrive at the critical curves in equations C.3.
In section 3.2, we saw how antisynchrony might become stable for large skewness in the absence of spike width effects. For finite spike width, this stability is still seen ( and in Figure 9a and in Figures 10 and 11) due to the dominant effect of the positive lobe of the PRC ( and in Figure 9c), but the region becomes sensitive (for type 1 PRCs) or disappears (for both types of PRCs) at large W/T. This observation is also verified numerically in Figures 9b and 9d.
For small skewness, the spike upstroke contributes a stabilizing eigenvalue component in both type 1 and type 2 PRCs, and the spike downstroke uniformly contributes a destabilizing component in both type 1 and type 2 PRCs because the sign of the PRC during the spiking of V(t−T/2) is not altered (see Figure 2B, case a, type I and type II). But as the skewness is increased, the spike upstrokes and downstrokes of V(t−T/2) begin to occur during the Z2(t) segment that is positive for type 1 and negative for type 2. Hence the role of spike upstrokes and downstrokes is reversed for type 2 PRCs (see Figure 2B, case e). Additionally, for type 1 PRCs, the boundary becomes a sensitive function of spike width because the segment Z2(t) is positive and the slope of the spike downstroke goes from infinity to a finite value quickly with W/T (see Figure 2B, case e, type I). Thus, the stable region bounded by and ( and in Figure 5a) shrinks drastically with increasing W/T (see Figure 10) for type 1 PRCs. As the skewness is increased, the contribution of the spike portion of the time course to the stability diminishes due to decreasing PRC levels, but the positive slope of the depolarizing phase that occurs during the negative PRC lobe contributes to the instability. Thus, eventually for type 2 PRCs, the stability region at large skewness deceases in size. Figure 10 also reveals that the antisynchrony that is stable due to large skewness totally disappears at very high frequency for type 1 PRCs, and is mostly confined to large negative B/C.
For the HH voltage model parameter of W/T=0.075, it can be verified that the stability diagram is similar to Figure 10b with little change in and , and the curves at the top left quadrant moving up slightly on the right and moving down slightly on the left. For the HH model's PRC parameters of A/T=0.567, and B/C=−0.5, the system lies in the bistable region above the curve and slightly to the right of . This is also verified numerically in Figure 1. The system is closer to the boundary of antisynchrony (), and thus altering the parameters in the HH model that solely reduces skewness is likely to cause the system to lose its stable antisynchronous state.
For the HH model's parameter a3/a2=0.2234, the stability diagram (for slightly bigger A/T) as a function of the product of spike with and frequency is almost identical to Figure 11d. At very low frequency, both synchrony and antisynchrony are stable in a bistable manner, and as the frequency is increased, only synchrony remains stable. We earlier remarked that for models such as LIF, a similar scenario happens, thus making the PRC skewness more critical in determining the stability states than its type. Many neuronal models are likely to have B/C ranging from −0.5 to 0.5, and in this range, as we can see from Figures 11c and 11d, PRC skewness plays an important role in stability. The PRC type plays a role in a qualitative manner, as, for example, the stability curves are not vertical to the B/C axis, but are sloped at some angle such that changing the type of the PRC may have a drastically different effect only at certain parameter points. A large negative B/C of course has different behavior from that of a large positive B/C.
The stability structure as a function of W/T and B/C is not altered significantly when the spike threshold is increased (see Figures 11e and 11f), but if the threshold is decreased to a very small value (see Figures 11a and 11b), the PRC type can cause qualitative changes in the stability structure. At such small spike thresholds, the depolarizing part (V2(t)) does not contribute much to the stability of either synchrony or antisynchrony. The spike upstrokes and downstrokes dominate the eigenvalue components. And when the skewness is large, the effect is more tangible: a portion of the curve is very close to B/C=0, and thus type 1 neurons at low frequencies show only synchrony, whereas type 2 neurons can show bistability. At moderately high frequencies, a portion of becomes very close to B/C=0 such that type 1 PRC neurons show neither synchrony nor antisynchrony, whereas type 2 PRC neurons show stable synchrony.
5.1. Relation to Previous Work
Synchrony among pulse-coupled neuronal oscillators received generic treatment earlier (Mirollo & Strogatz, 1990; Abramovich-Sivan & Akselrod, 1998; Goel & Ermentrout, 2002; Achuthan & Canavier, 2009) where the shape of the voltage was generalized by Mirollo and Strogatz, the slope of the PRC near zero phase was related to the stability of synchrony by Goel and Ermentrout and Achuthan and Canavier, and the shape of a type 1 PRC was parameterized to study the synchrony of nonidentical oscillators by Aromovich-Sivan and Akselrod. But synaptic or electrical coupling has always been difficult to generalize because of the inherent complexities of the models and the ensuing shapes of the voltages and PRCs. For electrical interactions, leaky integrate-and-fire (LIF) models became attractive in the past because for them, the voltage and PRC shapes can be derived analytically. Chow and Kopell (2000) and Lewis and Rinzel (2003) exhaustively studied the effect of oscillation frequency, and voltage shape on the phase-locked states using either the LIF model or modified versions of LIF model that incorporated spike and shape effects. Pfeuty et al. (2003) studied the effect of oscillation frequency and voltage shape parameter (that approximately corresponds to a3 in our model) in a quadratic LIF model. Lewis and Rinzel found that electrical coupling led to synchrony at all frequencies of the LIF. Integrate-and-fire neurons are of type 1 with zero spike width, and the PRCs computed for very small perturbations possess large skewness (see, e.g., Goel & Ermentrout, 2002). Our results for type 1 neurons in the absence of spike width and large skewness (see Figure 9) are consistent with their conclusions. Lewis and Rinzel also implemented spike effects by adding a positive term to the voltage profile every time a spike happens in the opposite neuron. This quantity emulates the missing spike profile in their model, and in our model, such a term is equivalent to considering spike upstroke effects. But during the spike upstroke, our PRC profile is zero; if it were nonzero, the synchrony region would be expanded to larger frequencies. In our model, synchrony could become unstable due to the spike downstroke effect that was not part of their model. Lewis and Rinzel also found that antisynchrony becomes unstable at high frequencies. Our bifurcation diagram at large skewness (see Figure 9) is in agreement with that observation, but we also find that had the skewness been smaller, then antisynchrony would have been stable at higher frequencies instead of at lower frequencies.
Chow and Kopell (2000) studied an LIF model that again is of type 1 but with a more complex structure to the voltage evolution that now incorporates not only spike amplitude but also spike width, slope of spike upstroke, and spike frequency. Such an incorporation was done as an added kernel to the voltage evolutions. Their kernel has the effect of altering both the voltage spike profile and the PRC at the same time as the kernel parameters are altered. In our formulation, we change them independently. But their general conclusion on the synchronous state is that it is stable at low frequencies and unstable at high frequencies. Our bifurcation diagrams for type 1 PRCs (see Figure 10) clearly map how such instability arises with increasing spike width or frequency. In Figure 11 the loss of synchrony is seen beyond a critical frequency. Chow and Kopell found the antisynchrony in their model at both high and low frequencies. Our results are at some variance with theirs. We find the antisynchrony stable either at low or high frequencies, but not in general at both the frequencies (see Figure 11). It is possible that the PRCs corresponding to Chow and Kopell's LIF model are acquiring distinctly different shapes from those of a pure LIF model because of the kernel affecting the voltage evolutions. At large skewness, for example, we do find within our model (see Figure 10) regimes of stable antisynchrony at low and high frequencies.
Pfeuty et al. (2003) studied a quadratic integrate-and-fire model that displays a type 1 PRC with more symmetry, albeit with jumps at the end, than that of an LIF model, but with small skewness. In particular, they reported that as the magnitude of the ratio of reset voltage (Vr) to the spike threshold voltage (VT) in increased stable synchrony is achieved, reducing this ratio leads to stable antisynchrony (with perhaps regimes of bistability with synchrony). The PRC does not have as much skewness as the LIF does, and hence falls under low skewness. We have not incorporated edge effects of the PRCs, if any, in our current study. But within our model, we see from Figure 8b that as the ratio a3/a2 is increased (i.e., the separation of reset and threshold levels increases and hence the magnitude of their ratio, Vr/VT), we find stable synchrony at large a3/a2, and stable antisynchrony at small a3/a2 and also at large frequencies.
Nomura et al. (2003) studied the Erisir et al. (1999) model and discovered that the interneurons are capable of displaying synchrony even at high frequencies. The model displays a type 1 PRC with small skewness, and the voltage time course has a finite spike width. This may appear counterintuitive because even in the Chow and Kopell model, synchrony became unstable at high frequencies, and in our results, it does become unstable for type 1 PRCs above a critical frequency (see Figure 11c) unless the PRC has large skewness (see Figure 11d). We wanted to investigate this puzzle and present the results in Figures 12a to 12f. The model (thin curves) does show stable synchrony (negative slopes of the growth functions) at high frequencies. But the PRC is no more of type 1 at high frequencies. The skewness is slightly increased (A/T from 0.17 to 0.27), but the type parameter (B/C) became negative (from 0.44 to −0.27) when the external steady current is increased to cause the frequency change from 64.3 Hz to 186.9 Hz. From Figure 11c we see that increasing the frequency while decreasing the type parameter keeps the system in stable synchronous state. We have also fitted this model's voltage (see Figures 12a to 12d) and PRC (see Figures 12b and 12e) with PWL functional forms (thick curves), and the resultant growth function (see Figures 12c and 12f) predicted the stability of the phase-locked states accurately at both these frequencies.
Mancilla et al. (2007) studied neocortical interneurons that exhibited type 1 PRCs at low frequencies. The PRCs retained their type 1 character even at 50 Hz. The resultant growth functions predicted stable synchrony and unstable antisynchrony. We have fitted piecewise linear curves to the average voltages and PRCs recorded from these neurons to determine our fit parameters (see Figures 12g to 12i). The spike width, skewness, and type parameters are all within the ranges of our model parameters. As the frequency is increased from 28 Hz to 50 Hz, the W/T factor is increased twice, the skewness is not altered much, but the type parameter B/C decreased from positive to 0. In both cases, the growth function computed from the PWL fits (see Figure 12i, lines) predicts the experimentally determined (see Figure 12i, points) stability of synchronous and antisynchronous states accurately. The reduction in the type parameter perhaps also indicates that increasing the frequency further might result in a type 2 PRC (with negative B/C) just as in the Erisir et al. model.
5.2. Effect of Spike Upstroke and Downstrokes
In type 1 PRC neurons, spike downstroke helps destabilize both synchronous and antisynchronous states (see Figure 2). In our PRC model, the segment at large phases is set to zero (Z5=0); thus, its contribution is not visible in the case of synchrony. But making this segment nonzero (but positive) causes both the synchrony and antisynchrony to enhance their ability to become stable due to the PRC being positive and the corresponding voltage segment also having a positive slope. The sign of the contribution of the upstroke segment () can also be inferred from Figure 2 by imagining the Z4(t) segment extending all the way to T. The sign would be the same as that of because the slope of the upstroke (V3(t)) has the same sign as that of the depolarizing phase (V2(t)). The spike upstroke helps stabilize the antisynchrony in all the cases considered (see Figure 2B).
In type 2 PRC neurons, spike upstroke always contributes to stabilizing the synchronous state (since the sign of is the same as , as explained above), but it can contribute to either stabilizing (cases a, b, j, k, l, and m in Figure 2B) or destabilizing (in the remaining cases) the antisynchronous state. That is, if the skewness is small such that A/T<(1−W/T)/2, then the upstroke contributes to stabilizing the antisynchronous state; otherwise, it contributes to destabilizing it. Spike downstroke contributes to stabilizing the synchronous state only if A>2W (i.e., in cases b and c); otherwise, its contribution toward stabilizing or destabilizing the synchrony depends on the sign of the PRC segment during the downstroke. Again in contrast, the spike downstroke contributes to destabilizing the antisynchronous state if the skewness is small such that A/T<(1−W/T)/2 (cases a, b, j, k, l, and m in Figure 2B). Similar behavior occurs in cases n and o. The downstroke contributes to stabilizing the antisynchrony at large skewness and small W/T such that A>T/2 and W<A/4 (cases d and e). In all other cases, the contribution of the spike downstroke toward stability of antisynchrony is dependent on the signs of the PRC segments during the upstroke.
5.3. Role of Spike Width and Spike Frequency
Type 1 PRC neurons are always in a stable synchronous state in the absence of the spike width effects for any level of PRC skewness (see Figure 5) because there is no factor that counters stable synchrony (see Figure 2A, case c, type I). But finite spike width, such that the product of the spike width and spike frequency (W/T) is large, can destabilize it (see in Figure 2A, case b, type I). For small or no skewness, instability occurs at smaller and smaller levels of W/T as the height of depolarizing voltage segment (a3) becomes smaller and smaller (, in Figure 8b). When the expression for is Taylor expanded, we get an approximate expression for synchrony when B=0 and A=0: Skewness helps increase the chances of synchrony in type 1 PRC neurons, competing against the influence of spike width. For finite skewness (A>0), sharp spike downstrokes are encountered for small spike widths. And when the spike downstroke does not extend beyond a factor that is proportional to the skewness A, synchrony becomes unstable (see Figure 10). But stability is retained for some values of B even when the skewness is smaller than 2W (curve in Figure 10) provided B/C is below the curve .
Type 1 PRC neurons can display antisynchrony due to either the dominant effect of spike upstroke or sharp PRC segments at large skewness. For small skewness, antisynchrony becomes stable only at sufficiently large frequencies (curves and in Figure 8). Again at A=0 and B=0, a polynomial can be fitted to giving an approximate condition for stable antisynchrony as Antisynchrony does not necessarily become stable immediately after the loss of stable synchrony. Other phase-locked states could mediate the in-between regime (see Figures 8d). As the skewness (A) is increased, the slope of the spike upstroke encounters smaller levels of the PRC, leading to its (upstroke's) diminishing effect on stabilizing the antisynchrony. Thus, the stable antisynchronous region that occurs at large frequencies (contributed to by the spike upstroke) decreases and disappears as the skewness is increased. As the skewness is increased, the slopes of the PRC segments to the left (Z3(t)) and to the right (Z4(t)) of the peak increase, leading to an increase of the negative eigenvalue components (, ). This results in stabilizing the antisynchrony even at low frequencies (W/T) (see Figure 9). But the slope of the left segment (Z3(t)) will decrease with increasing B, and thus its stabilizing effect diminishes. The stable antisynchronous region is then confined to smaller values of B/C at large skewness.
Type 2 PRC neurons, like those of type 1, can display both stable synchrony and stable antisynchrony. In addition to the spike width countering the synchrony as in type 1 PRC neurons, the PRC delay could also counter the synchrony (see Figure 2, cases b and c, type II). In the absence of spike effects (W=0), unlike type 1 PRC neurons, type 2 PRC neurons do not display synchrony for all levels of skewness. Synchrony fails at large skewness because of the increasingly destabilizing effect of the negative segments of the PRC ( and positive portion of in Figure 2, case c, type II). If B is large negative, these segments contribute more, and instability occurs at even smaller A. In particular, synchrony is stable (see equation 3.2 and Figure 5) only if So, for example, if maximum PRC delay is equal in magnitude to maximum PRC advancement, then synchrony fails if skewness is bigger than half of the spike period. Comparison of type I and type II plots in cases a, b, and d of Figure 2 reveals that some of the positive eigenvalue contribution that existed in type 1 neurons within the time span of spike downstroke is being removed in the type 2 neurons, and thus for small W/T, type 2 neurons are likely to remain synchronous in the presence of skewness. That indeed is the case, as can be seen in Figures 10a and 10b. However, as the spike width and frequency increase and if B is not sufficiently large negative, the gain in the form of the removed contributions referred to above can be countered by the effect of spike downstroke in combination with the positive segment of the PRC, leading to unstable synchrony at large W/T. Thus, in theory, we could have two boundaries for stable synchrony—one at small skewness and the other at large skewness. But for most of the negative B range, either the upper limit defined by and or the lower range defined by and constitutes the boundaries (see Figure 10).
Type 2 PRC neurons display antisynchrony, and its region in the parameter spaces can become fragmented into two, just as in the case of type 1. In the absence of spike width effects (see section 3.2), the spike downstroke effect inherent in the discontinuity of V(t−T/2) at t=T/2 counters the stability. The combined effect of the spike discontinuity (the sum of and , which is 2a3(AB+C(T−2A))/(T(T−A)) could contribute either a positive or a negative component to the eigenvalue. For A=0, this term is positive and helps destabilize the antisynchrony. For nonzero A, unless B is large negative (less than −C(T−2A)/A), the term continues to be positive and thus contributes to destabilizing the antisynchrony. This destabilizing effect is overcome by reducing the level of PRC at t=T/2, which can be achieved by increasing the skewness and pushing the PRC to the right. It requires a finite skewness () for stabilizing antisynchrony when W=0. But increasing the skewness to large levels has the opposite effect, because then the spike upstroke and downstroke fall in the negative lobe regime of the PRC, causing the upstroke to destabilize (due to its positive eigenvalue component) and the downstroke to stabilize in contrast to what happened at smaller A. Additionally, a considerable amount of the depolarizing phase could also fall in the negative PRC lobe regime. These factors together increase the chances of destabilizing the antisynchrony, and the instability occurs when (see Figure 5). We obtain stable antisynchrony when in the absence of spike width effects. Now consider the competing effects at a nonzero W (see Figure 2B, case b, type II) where the spike upstroke and downstroke assume less sharp slopes, but at the same time, the upstroke spreads toward higher levels of PRC and the downstroke toward lower levels of PRC. In effect, the effect of the upstroke (which contributes to instability) dominates, making it slightly harder to achieve stable antisynchrony. Thus, the stability will be lost as W/T is gradually increased (the stable region bounded by in Figure 11d), and also the level of A at which stability occurs (due to countering the spike effect by skewness) moves to larger values of A at finite W (, , , in Figure 10). Another mechanism could lead to stable antisynchrony at small A. When B is negative but not very large in magnitude, the reduction in the PRC level encountered by the spike upstroke may not be very significant (see Figure 2B, case a, type II). This results in a sizable stabilizing effect contributed by portions of the depolarizing phase and the spike upstroke to counter the destabilizing effect due to the spike downstroke and the negative PRC lobe. This results in a stable antisynchronous state occurring at small skewness levels. This effect lasts until the skewness sufficiently diminishes the PRC level during the spike upstroke (, , , in Figure 10).
5.4. Distinction Between Type 1 and Type 2 PRCs
This discussion makes it clear that both stable synchrony and stable antisynchrony occur in both type 1 and type 2 PRCs, but the amount of skewness, the type parameter, and the spike width or frequency control where these phase-locked states are located. For example, large positive B at large skewness and large negative B at small skewness hold stable synchrony in large parameter areas for some spike frequencies (see Figure 10c). Similarly, small skewness with large positive B, and large skewness with large negative B hold stable antisynchronous solutions (see Figures 10a to 10d).
But at the level of neuron pairs, any critical line that starts at B=0 and runs perpendicular to the B-axis would constitute a clear distinction between type 1 and type 2 PRCs in terms of their ability to display stable synchrony and antisynchrony. Such critical lines can be found ( and ) when the voltage time course has zero slope (when a3/a2=0, which is closer to the cases depicted in Figures 11a and 11b), where brief parameter ranges along W/T exist that clearly distinguish type 1 and type 2 PRCs. Consider the eigenvalue corresponding to case b for synchrony (). Examining the components (in appendix D), we find that only one component () becomes nonzero when a3=0, and it is equal to (curve ), which clearly changes its sign, causing a change of stability of synchrony when B is increased from negative to positive. This eigenvalue regime is defined (see the definition of Region b in appendix A) by the range . Rewriting this range, we conclude that when a3=0 and the frequency is high enough such that the spike width is bigger than a quarter of the skewness but less than half of skewness (), type 2 PRC neurons synchronize and type 1 PRC neurons lose synchrony.
Now consider the eigenvalue corresponding to case e for antisynchrony (). Examining its components (in appendix E) by setting a3=0 reveals that there are only two nonzero components ( and , corresponding to the spike upstrokes and downstrokes), which add up to . This clearly changes its sign as B is moved across 0 and is negative (leading to stability) for B<0 and positive (leading to instability) for B>0. The eigenvalue regime (see the definition of Region e in appendix B) is . Rewriting this range, we conclude that when a3=0, and when the spike frequency is small enough such that , antisynchrony is stable for type 1 PRCs and unstable for type 2 PRCs.
There is also a range of skewness for certain frequencies that comes close to distinguishing between type 1 and type 2 PRCs. The boundary of antisynchrony at large skewness for type 1 PRCs is sensitive to changes in frequency ( in Figure 10). For low frequencies around W/T=0.05, is positive and is near B=0. This clearly acts like a threshold for making the antisynchrony unstable.
In the specific cases mentioned, the parameter ranges near crossings of the critical boundaries across B=0 become regimes where crossing from type 1 to type 2 brings in qualitative change in the behavior of coupled neurons. Except in these cases, the stable boundary lines cross the B=0 axis with a nonzero angle, and thus in general, the PRC skewness and the voltage profile must also be considered in characterizing the PRC behaviors.
The main question we address in this study is when electrically coupled neurons synchronize under the assumption of weak coupling. Conventionally such a question is addressed by investigating the dynamics of the internal gating currents and then relating the mechanism to the emerging collective behavior of the coupled networks. Such investigative methods include drug application to block or enhance the activation of a selected ion channel type and application of steady current to cause a change in the neuron's oscillating frequency. These methods, in addition to causing the desired effect on physiological behavior of the neurons, also alter a number of other observables of interest such as spike width, spike frequency, spike height, PRC type, and PRC skewness. Although these observables are interrelated, insights into the effect of each of them on the network behavior could be detailed only if there was a mechanism to parameterize each of them. Such a task is difficult to carry out experimentally. But we can take advantage of the theory of weakly coupled oscillators to address this in a comprehensive manner.
Earlier studies that aimed at such generalizations largely focused on integrate-and-fire models (Chow & Kopell, 2000; Lewis & Rinzel, 2003; Pfeuty et al., 2003; Lewis & Skinner, 2012; Hansel, Mato, & Pfeuty, 2012) while principally working with fixed relationships between voltage time course and the PRC. In this study, we have taken the liberty of making these parameters independent (i.e., making the PRC and the voltage independent) and study the emerging collective states based on these shapes. This task would have been difficult but for the weakly coupled oscillator theory that is already being used widely in neuroscience (Van Vreeswijk et al., 1994; Hansel et al., 1995; Ermentrout, 1996; Netoff et al., 2005; Preyer & Butera, 2005; Ermentrout, Galán, & Urban, 2007; Tateno & Robinson, 2007; Cui, Carmen, & Butera Robert, 2009; Smeal, Ermentrout, & White, 2010). However, parameterizing the shapes independently would inevitably produce some parameter regimes that may be of less physiological interest, but at the same time it exhaustively covers all the possible and physiologically relevant voltage and PRC shapes within our model assumptions. The predictions would now become more general based on these two functions. In the current approach, only one independent time parameter (A/T, the PRC skewness) and one independent type parameter (B/C) define the PRC. The voltage spike train is again defined by just one independent time parameter (W/T, the spike width parameter), along with the spike peak, maximum spike hyperpolarization, and spike threshold level. But the latter three parameters collectively result in one truly independent parameter, a3/a2, that affects the stability of synchrony and antisynchrony. The stability is studied by segmenting the (W, A) parameter spaces (see Figure 2) such that each regime contains segments of the voltage and the PRC with constant slopes. This minimalistic formulation is possible because we have employed piecewise linear functions to model the voltage and the PRC shapes. Note, however, that the functional form of these shapes is less important than the ability to parameterize the shapes. We examined exhaustively whether the conventional categorization of PRCs into type 1 and type 2 had a bearing on their role in synchrony. More specifically, we addressed how the amplitude of the negative lobe in the PRC affects synchrony. We also investigated the role of PRC shape skewness as the maximum phase advancement of the PRC is systematically moved to longer phases. The effect of voltage depolarization and spike frequency was also investigated.
All of the analytical results relating the shape parameters are given in sections 3 and 4 and are summarized for zero spike width in Figure 5 and nonzero spike width in Figures 10 and 11. These results explain and put in perspective the previous results obtained using leaky integrate-and-fire (Chow & Kopell, 2000; Lewis & Rinzel, 2000), quadratic integrate-and-fire (Pfeuty et al., 2003), and Hodgkin-Huxley model equations (Nomura et al., 2003). The leaky integrate-and-fire model exhibits large PRC skewness, making a two-neuron LIF network exhibit antisynchrony at a small spike width and frequency (see Figure 9). The coupled network of oscillating HH model neurons exhibits bistability (see Figures 10b and 1) again because of the large skewness. Altering the PRC of the HH model such that the maximum PRC delay occurs at earlier phases would eliminate the antisynchrony. We also found (see Figure 11) that if the PRC skewness is small or the spike height is large (a3/a2 small), the synchronous state could be preserved even at high frequencies or large spike widths provided the type parameter B becomes sufficiently negative. This does occur in some models (Nomura et al., 2003).
In simulations, it may be possible to devise parameters or parameter combinations that could lead to independent control of the PRC and voltage shape parameters. But such a control may be more difficult in experimental preparations. For example, increasing a steady applied current or incrementally controlling the maximum ionic conductance of a selected channel could cause multiple changes in the PRC and voltage shapes, traversing some curvilinear path in our parameter spaces. But a simple assessment of the PRC and the voltages to acquire their shape parameters could be made. Our results provide guidance to the network behavior when such shape parameters are determined. The results when spike width is zero (see Figure 5) are the simplest to interpret. For PRCs whose normalized type parameter is bigger than −2(1−A/T) (curve in Figure 5), the coupled network displays stable synchrony. That is all type 1 PRCs lead to synchrony. Type 1 PRCs whose normalized skewness exceeds some threshold (see equations 3.15 to 3.18) lead to antisynchrony (curves and in Figure 5). Type 2 PRCs may cause loss of synchrony if the maximum PRC delay is large with sufficient skewness (curve in Figure 5). Large (bigger than ) but not very large (less than ) skewness in the type 2 PRCs causes stable antisynchrony. At zero spike width, the slope of the voltage depolarization does not determine the stability of synchrony or antisynchrony.
For finite spike width, the stability boundaries (see Figures 10 and 11) are more complicated. As discussed in section 5, the boundary of antisynchronous state, particularly when the PRC is of type 1, depends sensitively on spike width (curve in Figure 10). For sufficiently large spike width or at sufficiently high frequency, type 1 and type 2 PRCs can lead to an unstable synchrony, but the threshold of instability for the type 2 PRCs is higher than that for the type 1 PRCs. At high frequencies or large spike width, stable antisynchrony may disappear if the PRC has large skewness; instead, it may appear at small PRC skewness.
The limitations of our model come in the form of limitations of the voltage and PRC profiles and the basis of the weak coupling among oscillators. The voltage time course with three piecewise linear profile described HH time courses (see Figure 1), Erisir et al. (1999; see Figures 12a and d), and even some experimental recordings (see Figure 12g). The piecewise linear formulation of the PRC profile also described satisfactorily the PRCs of Hodgkin-Huxley (see Figure 1) and the experimental recordings of Mancilla et al. (2007; see Figure 12h), but some deviations are prominent between the piecewise curves and Erisir et al.’s model (see Figures 12b and 12d). These deviations are at longer input phases. Though these deviations had no effect in predicting the stability of the phase-locked states, it is also a reflection of the limitations of our five-segment piecewise linear PRC model in being unable to capture perhaps a large class of other PRCs. This model is also obviously unsuitable for PRCs with more than one peak (see some examples reported in Perez Velazquez et al., 2007; Devlin & Kay, 2001). In fact the profile of the PRC that we chose that has zeros at the ends may not be suitable for all experimental data. For example, a polynomial fit to the PRC that would have finite values at the ends may be more appropriate for some experimental recordings (Netoff et al., 2005; Tateno & Robinson, 2007). Also, because the coupling by electrical interactions uses all the voltage profile, deviations from the three piece-wise linear profile can have a significant effect on the stability of the phase-locked states, particularly when those deviations impart slopes with opposite signs. Finally, the interactions are assumed to be weak, and thus the results reported here cannot explain effects due to strong coupling. We have also not considered the role of heterogeneities and noise (for such studies in LIF models, see Ostojic et al., 2009) on phase-locked states. In fact the synchronous state considered here is a perfect synchronous state where the spikes align with no phase difference. The regions of unstable synchrony bordering the curves of stability contain phase-locked states that are near synchronous.
Appendix A: Computing the Eigenvalue for the Stability of Synchrony
We can use the profiles of Z(t) and V(t) to compute using equation 2.5. Since both Z(t) and V(t) are piecewise linear, the resultant will have piecewise contributions. As a typical case that we must deal with in computing the eigenvalue, consider the figure panel marked case c, type II in Figure 2A that corresponds to the arrangement of Z(t) and V(t) for models such as the HH model. The time point at which the spike downstroke ends (at t=2W) is less than the point at which the PRC becomes nonzero (t=A/2), and hence we obtain the condition, . The line A=4W is depicted in the (W, A) plot of Figure 2A, and the condition is represented by region c. The parameter point for the HH model itself is marked by a square symbol. There are six regimes that are separated by the change of slopes of these curves. Each regime contributes a portion (, , , , , and ) to the eigenvalue, and together all six regimes will combine to determine the stability of the synchronous state. In this case, we can directly see that , , and are 0 because Z(t) is zero in those regimes. The other eigenvalue components will have either positive or negative values depending on the product and its integral in those regimes. We shaded the regimes or parts of regimes that contribute to the positive value of the eigenvalue component. In this figure, we see that is completely positive because Z(t) (i.e., Z2(t) in particular) is negative in this regime and the slope of the voltage is positive. Hence the product is negative, resulting in a negative value of the integral over this regime. Since there is a minus sign already outside the integral in equation 2.5, the eigenvalue component becomes positive; hence, the region is shaded. In the adjacent regime that contributes to , the product changes its sign when Z(t) goes from negative to positive, and thus the positive eigenvalue portion that is marked by negative Z(t) is shaded. The sign of the component itself depends on the magnitude of the maximum delay B, the maximum advancement C, and the skewness A. These two components together determine whether the synchronous state can become unstable because the other components and are completely negative since the product is positive in those regimes. But if the PRC becomes type 1—in Figure 2a, case c, type I), all the eigenvalue components become negative, and thus the synchronous state cannot be destabilized.
As the skewness becomes smaller than 4W such that the maximum depolarizing time point (t=2W) is between the point at which the PRC becomes nonzero (t=A/2) and maximum delay (t=A), we obtain the condition (region b) in (W, A) plane of Figure 2A. The corresponding arrangements of Z(t) and V(t) are shown for type 1 and type 2 PRCs in case b of Figure 2A. When the product is negative, the resultant portion of the eigenvalue becomes positive, contributing to the instability of the synchronous state. Such regions along the time axis are shaded. In this case, both type 1 and type 2 can potentially cause the synchronous state to become unstable. Reducing the parameter A further such that the maximum delay (t=A) in the PRC occurs within the spike downstroke (t=2W) but the downstroke ends before the maximum phase advancement (t=A/2+T/2), we arrive at the condition (region a in (W, A) plane of Figure 2A). This condition occurs when the frequency is high. The corresponding arrangements of Z(t) and V(t) for type 2 and type 1 are shown, respectively, in case a, type II, and case a, type I, in Figure 2A. The regions that contribute to positive eigenvalue components are shaded. Finally, when the frequency of oscillation is very high such that the spike downstroke (t=2W) ends after the PRC acquires its maximum advancement (t=A/2+T/2), we arrive at condition A<4W−T (region d in the (W, A) plane of Figure 2A). An example of the arrangement of Z(t) and V(t) is displayed for type 2 and type 1 PRCs in case d, type II, and case d, type I, of Figure 2A. There are now larger regions shaded that mark the regions contributing to positive eigenvalue, which destabilizes synchronous state. We summarize the four parameter regions discussed above as follows:
- Region a:
Intermediate skewness: :
- Region b:
Large skewness: :
- Region c:
Very large skewness: :
- Region d:
Small skewness: :
Appendix B: Computing the Eigenvalue for the Stability of Antisynchrony
The eigenvalue for the antisynchronous state is found using equation 2.6, but this case turns out to be more complex than that for . comprises piecewise contributions from Z(t) and V(t−T/2). Depending on the value of the skewness (A) in relation to the spike width parameter W, the relative arrangement of Z(t) and V(t−T/2) changes while keeping the spike peak position of the voltage time course fixed at T/2. Four such arrangements are displayed in Figure 2B. At any given such arrangement of Z(t) and V(t−T/2), we will have to evaluate eight integrals corresponding to the eight different time segments during which the slopes of Z(t) and V(t−T/2) are constant. We refer to these integrals as eigenvalue components.
At A=0, for example, the spike peak and the positive peak of the PRC coincide, and for small positive A, the regions marked a and b in the (W, A) plane of Figure 2B arise. Examples of the arrangement of the PRC and V(t) in these two regions are shown, respectively, in cases a and b. Consider case a where the positive peak position () of the PRC is slightly shifted toward longer phases, but the spike width is still big enough that the peak occurs during the downstroke of the spike and before the spike minimum (at ). That is, we have the relation A<4W, and if A is increased beyond 4W, another set of regions along the time axis contributes to the eigenvalue. But W is small enough that the spike minimum (at ) occurred within the Z4(t) branch and before Z5(t) began (at ). That is, we must have . And A is also small enough that the position of PRC minimum (at A) occurs before the spike upstroke initiates (at ). Hence, we also have the condition that . If these conditions are violated, we will have to rewrite the conditions for the validity of the corresponding arrangement of Z(t) and V(t−T/2). These three conditions are valid in region a in the (W, A) space of Figure 2B, and the region is bounded by A=4W, W=T/5, and A=(T−W)/2, which are also drawn in the figure. Case a in the type II and type I panels also shows the regimes during which the slopes of neither Z(t) nor V(t−T/2) change. These regimes from the left (starting at t=0) contribute to the eigenvalue components . Since the slopes of V(t−T/2) and Z(t) are constant, the sign of the eigenvalue component can be directly visualized from the panels. When the product is negative either because of the negative voltage slope or because of the negative PRC segment, then it contributes to a positive eigenvalue (due to the extra negative sign outside the integral in equation 2.6). In type 2, the negative part of the PRC contributes to positive eigenvalues ( and part of are positive). In both type 1 and type 2 of case a, the spike downstroke contributes to a positive eigenvalue: and both become positive. The spike downstroke, even in the form of a voltage discontinuity when the spike width is zero, contributes significant positiveness to the eigenvalue due to the sharp voltage gradient. We will see that the antisynchronous state is unstable for small W even when A=0, and it becomes stable only for sufficiently long W that the downstroke contribution diminishes due to the decreased slope.
Starting from the conditions of region a, when A is increased slightly, we will arrive at three other possible arrangements of Z(t) and V(t−T/2). We will arrive at region b when the maximum PRC (at ) occurs after the spike minimum at , which results in condition A>4W while still obeying . Instead we could have the PRC minimum (at A) occur after the spike upstroke event (at ) but before the spike peak (at ), thus defining the bottom and top boundaries of region f: and . As a third possibility we could have both the above occurring at the same time: the PRC minimum (at A) occurring after the spike upstroke initiation (at ) but before the spike peak (at ) and the PRC maximum occurring after the spike minimum (at ). Together this gives the region of validity of this scenario at the intersection of the previous two regions and results in region c. In the example of case b, the regions along the time axis that contribute to a positive eigenvalue component on account of the product becoming negative are shaded in Figure 2B, case b.
The arrangements of Z(t) and V(t−T/2) in two (c and f) of the above four regions define the maximum value of A, which is smaller than on account of the PRC minimum (at A) occurring before the half-period (at T/2). The HH model is an example of its PRC minimum occurring above its half-period point (the filled square in region d) in the (W, A) panel of Figure 2B. But in general when A is increased further such that the PRC minimum occurs after the half-period point, the arrangement of the PRC and Z(t) in region c transforms into that of region d and region f into that of region g. And when A increased further such that the PRC minimum (A) occurs not only after the half-period point as shown in case d but after the spike minimum (at ), then we arrive at region e, for which case e shows an example. All of the seven regions occur when the spike width is less than one-fifth of the time period, . The arrangements of the PRC and V(t) in regions d and e are illustrated, respectively, in Figure 2B, cases d and e. The conditions for these seven regions are summarized below:
- Region a:
- Region b:
- Region c:
- Region d:
- Region e:
- Region f:
- Region g:
Regions a, f, and g border the maximum W that can occur in all seven regions by requiring , and other regions also obey this relationship. This condition arises because of the existence of the zero phase regime in the PRC at later phases, Z5(t). In all three regions, the spike minimum (at ) occurred before the Z5(t) segment began. When the spike width is increased such that the spike minimum occurred within the Z5(t) segment of the PRC (i.e., after and before t=T), then the arrangement of the PRC and V(t) in regions a, f, and g, respectively, become candidates for regions j, i, and h. The conditions for these three regions are given by:
- Region h:
- Region i:
- Region j:
The extremely large spike width with W above leads to some more interesting cases. Consider the parameter region when the spike upstroke (at ) occurs after the PRC minimum (at A) in segment Z3(t) as in regions a, b, and j: . If the spike width grows long enough such that it is now bigger than a quarter of the period, the spike minimum of V(t−T/2) will occur beyond the Z3(t) segment—in segment Z1(t) (resulting in region k), segment Z2(t) (resulting in region l), or segment Z3(t) itself (resulting in region m). If the spike upstroke occurred in segment Z2(t) instead of Z3(t) because of an elevated A that moves the PRC to longer phases, then for long spike widths, there will be two cases corresponding to the discussion above: one when the spike minimum (now at ) occurs in segment Z1(t) (resulting in region n) and the other when the spike minimum occurs in segment Z2(t) (resulting in region o). And finally when the skewness is large enough that the spike peak occurs before the PRC minimum but within segment Z2(t), then for large enough spike widths, the minimum of V(t−T/2) can occur in either the PRC segment of Z1(t) leading to region p or the segment of Z2(t) itself leading to region q. All seven regions occur when the spike width is bigger than but still smaller than , as defined by equation 2.8. The conditions for the seven regions are given by
- Region k:
- Region l:
- Region m:
- Region n:
- Region o:
- Region p:
- Region q:
The spike downstroke or, when the spike width is zero, the discontinuity associated with the spike due to its sharp slope contributes maximally to the positive eigenvalue and can by itself have a larger magnitude over the sum of all the other negative eigenvalue components. So even when there is no spike width (W=0), the antisynchronous state can become unstable due to the sharp voltage drop at t=T/2. When the spike width increases, its effect in destabilizing the antisynchrony decreases, and the antisynchronous state has a greater chance of becoming stable.
Appendix C: Boundaries of Antisynchrony When W>0
Appendix D: Eigenvalue Components
D.1. Region a
D.2. Region b
D.3. Region c
D.4. Region d
Appendix E: Eigenvalue Components
E.1. Region a
E.2. Region b
E.3. Region c
E.4. Region d
E.5. Region e
E.6. Region f
The first and last components are zeros: and . The other components are identical to previously computed expressions in either region a or region c: