Abstract

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.

1.  Introduction

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

Coupled oscillatory model neurons that are intrinsically nonlinear in nature (such as those modeled by Hodgkin-Huxley equations) can be described by coupled phase evolution equations applying the theory of weakly coupled oscillators. The main idea behind such a reduction is that one would be able to study the collective behavior of coupled neurons with information that is easily obtainable in experiments. Instead of the full intrinsic dynamics that are necessary to model individual neurons, the phase-coupled model requires only the knowledge of phase response curves, the voltage time course, and the coupling mechanism, all of which are relatively easy to obtain. If and represent phases (i.e., times elapsed since their last respective spike times) of the two coupled neurons that are oscillating with nearly identical periods T1 and T2, their collective behavior is described by the two-phase evolution equations (Winfree, 1967; Neu, 1979b, 1979a; Kuramoto, 1984; Ermentrout & Kopell, 1986, 1991; Hoppensteadt & Izhikevich, 1997; Brown et al., 2004):
formula
2.1
where is a small (because of the assumption of weak-coupling) constant parameter, and are the mutual interaction functions that are essentially the instantaneous frequency increments. We assume that the oscillators are identical, and thus T1=T2=T, and , where is given by a convolution of the PRC, Z(t), that is characteristic of each of the oscillators, and the coupling function, which in our case of electrical coupling is simply the difference of the voltages [V(t)] of the two neurons:
formula
2.2
If is the interaction function of the first oscillator, the second oscillator's interaction is defined by . The voltage with phase advancement [] is from the coupled oscillator and is also the factor that will make the computation of the stability a complex task even within our assumption of piecewise linear functions. The dynamics of the two identical coupled oscillators can be reduced further to the study of a single equation by writing an equation for the phase difference as
formula
2.3
where
formula
2.4
We will call the function the growth function because it essentially quantifies the growth of the phase difference, or the rate of divergence of the phases. If there is no divergence, that is, when at some , those phases represent equilibria of the original phase equations. These equilibria can be synchronous states (), antisynchronous (), or any other phase-locked states. We can clearly see from equation 2.4 that becomes zero at and T/2, and thus synchronous and antisynchronous states always exist for this model. We are concerned in this study with these two states. However, to be of any practical utility, these states must be shown to be linearly stable. An equilibrium is stable if any perturbation imparted to the system subsides in time, that is, the derivative of the growth function at the equilibrium must be negative or the eigenvalue , where , must be negative. We represent this eigenvalue by if , and if . We can directly express these eigenvalues from equation 2.3 and 2.4 as
formula
2.5
formula
2.6
In this study we compute and for various shapes of Z(t) and V(t), and map regions of stability of both synchronous and antisynchronous states. But first we formulate V(t) and Z(t).

2.1.  Model for the Voltage Time Course

The voltage trace V(t) is periodic in time with period T (>0) and is formulated by three piecewise linear curves (see Figure 1) that use three parameters (Vp, Vm, and Vth) defining different voltage levels, and another parameter that is closely related to the width of the spike that is asymmetric. This formulation is inspired by an empirical observation of the time course of the Hodgkin-Huxley (HH) model neuron (depicted as thin lines in Figure 1 for an applied current of A/cm2). The three linear curves correspond to the spike upstroke, downstroke, and depolarization regimes. The number of free parameters is kept to a minimum, and they are Vp, the peak value of the spike; Vm, the maximum hyperpolarization reached by the downstroke of the action potential; Vth, the spike threshold; and finally the parameter W. These parameters are now freely variable and no longer correspond to the HH model alone. The voltage profile is thus given by
formula
2.7
By observing the spike profile and the time course, we arrive at some simple conditions on the parameters defining V(t), as described below. The spike width and the refractory period together are accounted for by the rise time (W/2) and the fall time (2W), and we insist that the spike period (T) be not smaller than this so that there is a finite time to recover before the next spike onset: , or
formula
2.8
The actual width of the spike could be considered either 5W/2 or less, but for simplicity, we term W in this study as the spike width. The spike peak Vp is bigger than the Vth, and hence
formula
2.9
The spike peak Vp can be either positive or negative, but its magnitude is assumed to be not smaller than that of Vm, hence
formula
2.10
The maximum hyperpolarization Vm is assumed to be negative, and the magnitude of Vth is smaller than that of Vm, such that
formula
2.11
The shape of the voltage trace also imposes the condition
formula
2.12
and
formula
2.13
Of all the parameters that define the voltage shape, we will see that only W/T and a3/a2 are the truly independent parameters that affect the stability boundaries of synchrony and antisynchrony.
Figure 1:

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

Figure 1:

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

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.

This formulation, as in the previous section, is again motivated by an empirical observation of the PRC shape of the HH model (thin lines in Figure 1). Thus, the piecewise linear formulation incorporates some general features of the HH model, but since the shape parameters are freely variable, our predictions based on the piecewise linear model are applicable to a much larger set of models than those based on the HH model alone. The profile of Z(t) is given by
formula
2.14
and is illustrated in Figure 1 (thick lines). The parameter B is the type parameter that can convert the PRC to become completely non-negative () and thus a type 1, or partially negative (B<0), and thus a type 2. We term A the skewness parameter because the nearness of the maximum phase delay position to the center of the period is controlled by A. By this definition, A=0 signifies a more symmetric PRC that has the maximum phase advancement at half-period. Since the maximum advancement occurs at and should be less than above which the PRC is zero, this leads to the following condition:
formula
2.15
The parameter C(>0) is the maximum phase advancement of the PRC. We will see that B/C and A/T are the only independent PRC parameters that affect the stability boundaries of synchrony and antisynchrony.

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

In the synchronous state, the phases of the two coupled oscillators are identical: the spike time difference of the two neurons is zero, and in the antisynchronous state, the phase difference is equivalent to half the oscillation period. The stability of these two states is determined, respectively, by and defined in equations 2.5 and 2.6. In computing and from these formulas, it is convenient to segment the parameter space of (W, A) depending on the relation between the V and Z segments. We have detailed this procedure in appendix  A. To help us compute the stability of synchronous state, the (W, A) parameter regime is split into four regions: a, b, c, and d. These regions and the relative position of the V and Z in each of these regions are displayed in Figure 2A. The eigenvalue in each region is written as
formula
2.16
The expressions for each of these components are listed in appendix  D. To help us compute the stability of antisynchronous state (i.e., find from the integral in equation 2.6), we time-shift the voltage by a half-period, and identify the parameter regimes where slopes are constant. This results in 17 different regions (a, b, c, , q), and the eigenvalue in each region is written as
formula
2.17
Figure 2:

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

Figure 2:

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

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.

Figure 3:

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

Figure 3:

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

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.

Figure 4:

When the spike width is zero, a large skewness in type 2 PRCs can destabilize a synchronous state, still cause bistability between synchronous and antisynchronous states, and destabilize both synchronous and antisynchronous states. (a) Voltage time course from equation 2.7 at W/T=0. (b) PRC profiles from equation 2.14 when B/C=−0.5. The eigenvalue and the different components that make up the eigenvalue, which determines the stability of the synchronous (c) and antisynchronous (d) states. The bifurcation diagram (e) shows stable (solid curves) and unstable (dashed curves) phase-locked states as the skewness is increased. Different stable regimes are marked pictorially in panel f.

Figure 4:

When the spike width is zero, a large skewness in type 2 PRCs can destabilize a synchronous state, still cause bistability between synchronous and antisynchronous states, and destabilize both synchronous and antisynchronous states. (a) Voltage time course from equation 2.7 at W/T=0. (b) PRC profiles from equation 2.14 when B/C=−0.5. The eigenvalue and the different components that make up the eigenvalue, which determines the stability of the synchronous (c) and antisynchronous (d) states. The bifurcation diagram (e) shows stable (solid curves) and unstable (dashed curves) phase-locked states as the skewness is increased. Different stable regimes are marked pictorially in panel f.

Figure 5:

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

Figure 5:

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

3.1.  Synchrony

Neurons with type 1 PRCs that have finite zero segments at either end always display synchrony as there are no eigenvalue components that can make the synchrony unstable (see Figure 2A, case c, type I). But if the PRC is of type 2, then synchrony can become unstable due to the negative PRC regime (see Figure 2A, case c, type II). Since W=0, the only eigenvalue components that are nonzero are , , which may become positive for large B, and (see Figure 4c). Combining these components, we get the eigenvalue We see that (see Figure 3c), and hence the synchrony is stable at B=0. A critical state of stability will be reached when becomes 0. By setting the eigenvalue to zero, we obtain . And on this critical curve, . This critical value acts as a lower boundary for stability. Increasing B across causes the eigenvalue go through 0 with a negative slope, and hence the eigenvalue is positive for and negative for and is zero at . Thus, the region of stable synchrony is given by
formula
3.1
where we have also normalized B with C. This curve is negative for A between 0 and T. Consequently, all type 1 PRCs that have zero segments at the edges display stable synchrony (see Figure 3f), but type 2 PRCs (B<0) do not necessarily lead to unstable synchrony. The maximum delay of the PRC (B) must be sufficiently long to counter the skewness in order to destabilize the synchrony. But PRCs with large skewness are closer to the unstable boundary than those with small skewness.
The above critical condition for stability can be inverted to obtain an expression in terms of the skewness. Thus, the stable synchrony is obtained when
formula
3.2
The above condition, shown in Figure 4f, predicts the value of the upper boundary on skewness for synchrony to be stable. For the HH model discussed earlier, B/C=−0.5, and hence the upper limit of skewness for synchrony is A/T=0.75. Numerical bifurcation diagrams were computed using functional forms of V(t) and Z(t) and are illustrated for type 1 (see Figure 3e) and type 2 (see Figure 4e), which verify the above analytical prediction.

3.2.  Antisynchrony

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.

First consider the case of . The eigenvalue components consist of , , , , and the integral across the voltage discontinuity at that can also be obtained by combining and in the limit of zero spike width. Combining these components and using the formula in equation 2.17, we write the eigenvalue as We directly see that for B=0, the antisynchronous state is unstable for small skewness and becomes stable for large skewness because
formula
Thus, the criterion for stability when B=0 is
formula
3.3
and is illustrated in Figure 3f. As B is increased or decreased from 0, antisynchrony may become stable. The sign of could change across a critical curve (), which is obtained by solving the equation ,
formula
3.4
where C is used to normalize B. On this critical curve, the transition of the eigenvalue is given by
formula
Thus, for very small skewness levels (), the eigenvalue becomes negative when B is increased past . Thus, the region of stable antisynchrony is given by
formula
3.5
where is given by equation 3.4. So when there is no skewness in the PRC (A=0), the type parameter must be two times larger than the maximum PRC advancement (B>2C) for the antisynchrony to become stable. For slightly larger A, the stability boundary increases quadratically, as seen in equation 3.4. This curve, , appears in the top left corner of Figure 5a.
For , the stable synchronous region lies below B=0, and crossing from below the eigenvalue becomes positive. Hence the stability region is defined by
formula
3.6
where is given by equation 3.4. This curve is shown in Figure 5a.
In the region , the antisynchronous state is already stable when B=0, and across the critical curve , it loses stability. Hence, the curve lies above B=0. Thus, the stability region is again given by
formula
3.7
where is given by equation 3.4. The curve and the region represented by the above relation is shown in Figure 5a.
Next consider the case of where the eigenvalue components (see Figure 2B, case e) are given by , , , , and the contribution from the voltage discontinuity that can in turn be computed by taking the limit of spike width going to zero from and . Combining all the eigenvalue components and using the formula in equation 2.17, we get the total eigenvalue as This eigenvalue is negative at B=0 in this regime because . Thus, the antisynchronous state is already stable. A critical curve that may lie above or below the B=0 level is obtained by solving for B from the equation , which gives
formula
3.8
This clearly lies above B=0 when , and below B=0 level when . The transition of the eigenvalue is determined by
formula
Thus, when , the stable antisynchronous state becomes unstable when B is increased above . The stable region for antisynchrony is given by
formula
3.9
where is given by equation 3.8. This region is shown in Figure 5a.
When , because the eigenvalue transitions to negative value across , the stable region lies above the critical curve, and thus the stable region for antisynchrony is given by
formula
3.10
where is again given by equation 3.8. This region is shown bounded by the curve in Figure 5a. Thus, type 2 PRCs with very large skewness can lose both synchrony and antisynchrony.
The critical curves can also be written in terms of A, and the regions can be inferred by solving the equations and . We list the critical curves and the regions below but infer the regions by simply matching them with those already derived above instead of repeating the analysis. For a specific value of B/C=−0.5, the eigenvalue components and the total eigenvalue are shown in Figure 4d and the stable antisynchrony region derived from the eigenvalue is marked in Figure 4f. The computed numerical bifurcation diagram, Figure 4e, verifies this region. From the equations for , critical curves are obtained as
formula
3.11
formula
3.12
And from the equation for , the following critical curves are obtained:
formula
3.13
formula
3.14
We list next the critical regions in terms of A that can be verified with those already given above by numerically plotting them. The antisynchrony is stable in the following regions:
formula
3.15
formula
3.16
formula
3.17
formula
3.18
For B/C=−0.5 as in Figure 4, the last condition gives us the stable range of skewness as 0.347<A/T<0.883. This region is depicted in Figure 4f.

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.

Figure 6:

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

Figure 6:

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

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.

Figure 7:

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

Figure 7:

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

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.

Figure 8:

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

Figure 8:

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

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.

Figure 9:

Large skewness (A/T=0.6 here) can make the antisynchronous state accessible for a range of W/T values starting at zero. Such a regime is thinner in type 1 than in type 2 PRCs. (a) Stability regions in (W/T, a3/a2) space for B/C=0.2, 0, −0.25. The curves and are the boundaries for synchrony. The curves and are the boundaries for the antisynchrony. (b) One-parameter bifurcation diagram showing stable (solid) and unstable (dashed) states as a function of W/T at a3/a2=0.22 that shows bistability of synchronous (phase-locked state at 0) and antisynchronous state (phase-locked state at T/2) at small W/T. (c) The eigenvalue components (thin curves) and the total eigenvalue (thick curve) that determine the stability of antisynchronous state as a function of W/T at a3/a2=0.22. (d) Growth functions at a3/a2=0.22 at a few values of W/T depicting how the antisynchronous state becomes unstable at large W/T.

Figure 9:

Large skewness (A/T=0.6 here) can make the antisynchronous state accessible for a range of W/T values starting at zero. Such a regime is thinner in type 1 than in type 2 PRCs. (a) Stability regions in (W/T, a3/a2) space for B/C=0.2, 0, −0.25. The curves and are the boundaries for synchrony. The curves and are the boundaries for the antisynchrony. (b) One-parameter bifurcation diagram showing stable (solid) and unstable (dashed) states as a function of W/T at a3/a2=0.22 that shows bistability of synchronous (phase-locked state at 0) and antisynchronous state (phase-locked state at T/2) at small W/T. (c) The eigenvalue components (thin curves) and the total eigenvalue (thick curve) that determine the stability of antisynchronous state as a function of W/T at a3/a2=0.22. (d) Growth functions at a3/a2=0.22 at a few values of W/T depicting how the antisynchronous state becomes unstable at large W/T.

4.1.  Synchrony

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.

Figure 10:

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

Figure 10:

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

Figure 11:

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

Figure 11:

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

4.1.1.  Very Large Skewness (A≥ 4W), Case c

This region extends all the way to A=0 but can be termed the case of very large skewness because the region extends beyond the maximum skewness in the other three cases. The Hodgkin-Huxley model (see Figure 1) falls in this regime. The arrangement of Z(t) and V(t) is shown for type 1 and type 2 in case c of Figure 2A. The spike profile does not affect the stability because it occurs during the zero phases of the PRC. Only three eigenvalue components are nonzero: , , and (computed in appendix  D). Combining these three and using the formula in equation 2.16, we get the eigenvalue as
formula
Utilizing the conditions in equations 2.15 and 2.8, we directly see that for . That is, if the parameter W/T is small enough or if the skewness is large enough such that 4W<A, then all type 1 PRCs lead to synchronous state and no instability occurs at larger A. In fact, synchrony is also stable for type 2 PRCs, and the condition for stable synchrony is obtained directly from the above equation when we insist that , which results in the following region:
formula
4.1
Whenever the spike width is less than one-quarter of the skewness, the above condition provides the lower boundary for synchrony, and there is no upper boundary. As the spike width is increased the regime of existence of the curve decreases, as can be seen from the panels in Figure 10. Since the skewness is limited by TW (see equation 2.15), as the spike width is increased, the skewness regime for the existence of falls in the forbidden region, and hence this curve will no longer define the boundary seen in Figure 10d. On the other hand, from the condition where exists, we see that as the skewness is increased, the span of W/T increases, as is also seen in the (W/T, B/C) plots in Figure 11. We can express the condition in terms of the skewness as well. From the eigenvalue equation, we see that it becomes negative if A<TW/2+BT/(2C). Given that A must be bigger than 4W and also smaller than TW, we write the sufficient condition for stable synchrony as
formula
4.2
This condition complements other conditions from the other three cases in defining the complete picture of stability. For the example shown in Figure 6, the parameter values are B/C=0, W/T=0.15, and thus the lower and upper limits in the above relation become 0.6 and 0.85. Since 0.6 is the lower limit of the present case, the lower boundary will have been complemented with the results from the other cases, but the upper limit on A/T is 0.85 (see Figure 6e). For the example shown in Figure 7, B/C=−0.5 and W/T=0.15, and thus the lower and upper limits in the above relation become 0.6 and 0.675. Thus, the upper boundary for synchrony is defined by A<0.675 (see Figure 7e).

4.1.2.  Large Skewness (2W≤ A> 4W), Case b

The arrangement of V(t) and Z(t) is shown for type 1 and type 2 PRCs, respectively, in case b, type I and case b, type II of Figure 2A. The nonzero eigenvalue components are , , , and and are listed in appendix  D. The critical curve is obtained by solving the equation that results in the following equation in terms of the normalized-type parameter:
formula
4.3
where , and the derivative of the eigenvalue on the critical boundary is
formula
By observing the behavior in the (W/T, A/T) plane (see Figure 10), we conclude that there are two segments to the curve , and they form boundaries of the synchronous state. These two segments fall on either side of the singularity, that is, between the segments. The right segment lies in the regime where B/C<−2+9W/T, and the left segment lies in the regime where B/C>e2; both segments reach a singularity when the denominator of the reaches 0, which is when the abscissa, A/T, asymptotes the following value:
formula
where , , , , and . Though the expression for is exact, it is difficult to see the dominant terms. Using the above asymptote, we get simple sufficient conditions for stable synchrony in this regime:
formula
4.4
formula
4.5
The asymptote can be seen in all the panels in Figure 10. For a3/a2=0.2234, using the boundaries predicted by the asymptotes for W/T=0.02, we obtain the stability boundaries as when B/C>0.58, and when B/C<−1.82. For W/T=0.05, the stability boundaries are obtained as when B/C>0.58 and when B/C<−1.55. For W/T=0.15, the stability boundaries are obtained as when B/C>0.6 and when B/C<−0.65. And finally for W/T=0.3, the asymptote results in a value of 0.7, which is also the boundary (1−W/T) of the existence of skewness, but the actual curve shows a small region of stable synchrony at this level of W/T. Thus, the asymptote analysis predicts the boundaries for small W/T very well, and for large W/T, the full expression for must be used.

4.1.3.  Intermediate Skewness (4W-T ≤ A < 2W), Case a

The arrangement of V(t) and Z(t) for this case is shown for type 1 and type 2 PRCs, respectively, in case a, type I and case a, type II of Figure 2A. The nonzero eigenvalue components are , , , and . The critical curve is obtained by solving the equation and is given in terms of B as
formula
4.6
where , , , , . The derivative of the eigenvalue at is
formula

4.1.4.  Small Skewness (0 ≤ A < 4W-T), Case d

And finally in this case, the arrangement of V(t) and Z(t) is shown for type 1 and type 2 PRCs, respectively, in case d, type I and case d, type II of Figure 2A. The nonzero eigenvalue components are , , , and . The critical curve that is obtained by solving is given in terms of B as
formula
4.7
where . The transition of the eigenvalue is given by .

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.

4.2.  Antisynchrony

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<(TW)/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(tT/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(tT/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(tT/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.  Discussion

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.

Figure 12:

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

Figure 12:

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

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(tT/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(TA)) 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.

6.  Conclusion

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<4WT (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: :

where is given in equation 2.16. The expressions for , , and are listed in Appendix  D for all four cases. The value corresponds to all or part of Z1(t), and corresponds to the contribution of Z5(t). Hence both of these are zero for all the cases. Consequently, the rising phase of the action potential has no effect on the stability of synchrony in our model. If we were to consider the PRC to be nonzero during the rising phase, that would contribute to a more stabilizing eigenvalue because the slope of the voltage is positive and PRC is always positive in that region. The effect of the falling phase of the action potential will have a diminishing effect as the skewness of the PRC increases. The parameter region of W above 2T/5 is forbidden by the condition in equation 2.8.

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(tT/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(tT/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(tT/2), we will have to evaluate eight integrals corresponding to the eight different time segments during which the slopes of Z(t) and V(tT/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(tT/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=(TW)/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(tT/2) change. These regimes from the left (starting at t=0) contribute to the eigenvalue components . Since the slopes of V(tT/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(tT/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(tT/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: 

    :

where is defined by equation 2.17. The eigenvalue components in equation 2.17 are computed using the same integration formula as in equation 2.6 with the limits adjusted according to regions a to g. The algebraic expressions for are listed in appendix  E. The values and are contributed by either part or all of the PRC segments Z1(t) and Z5(t), respectively, and thus they are zero. The other segments contribute either positively or negatively to the eigenvalue and decide the critical boundary that defines the stability region of the antisynchronous state.

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: 

    , :

where , x=h, i, j is the sum of the eight individual eigenvalue components in each region. The sum is as defined in equation 2.17, and the expressions for the individual segments are given in appendix  E.

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(tT/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(tT/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: 

    , :

where is the sum of the eight individual components as before as defined in equation 2.17. The expressions for the individual components are given in appendix  E. The region that is not displayed in the (W, A) plot in Figure 2B is forbidden due to the condition in equation 2.8.

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

The critical curves defining the antisynchronous state discussed in section 4.2 are listed below.
formula
C.1
The coefficients on the right-hand sides of the above expressions and similar expressions below are listed in appendix  F:
formula
C.2
formula
C.3

Appendix D:  Eigenvalue Components

D.1.  Region a

The following are the expressions for the six components of the eigenvalue that determines the stability of synchronous state under case a in appendix  A:
formula

D.2.  Region b

As in region a, the first and the last component become zero, and the fifth component is identical to that in region a: , , and . The other three components are:
formula

D.3.  Region c

The first two and the last component become zero: , , and . The fourth component is identical to that in region b, and the fifth component is identical to that in region a: , and . The other component is
formula

D.4.  Region d

The first and the last components become zero: , and . The second component is identical to that in region a: . The other three components are:
formula

Appendix E:  Eigenvalue Components

The expressions for the eight components (see equation 2.17) comprising the eigenvalue in each of the 17 parameter regimes in the (W, A) plane (under region a,, region q) for determining the stability of antisynchronous state described in appendix  B are listed below.

E.1.  Region a

formula

E.2.  Region b

Two of the components become zero, and three others are identical to those in region a: , and , , , and . The other three components are:
formula

E.3.  Region c

The components and . Three other components are identical to those in region b: , , and . The other three components are:
formula

E.4.  Region d

The first and last components are zeros: and . One component is identical to that in region c: . Two others are identical to those in region b: , and . The other three components are:
formula

E.5.  Region e

The first and last components are zeros: and . The third and seventh components are identical to, respectively, those in regions d and b: and . The other four components are:
formula

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: