Abstract

The role of the phase response curve (PRC) shape on the synchrony of synaptically coupled oscillating neurons is examined. If the PRC is independent of the phase, because of the synaptic form of the coupling, synchrony is found to be stable for both excitatory and inhibitory coupling at all rates, whereas the antisynchrony becomes stable at low rates. A faster synaptic rise helps extend the stability of antisynchrony to higher rates. If the PRC is not constant but has a profile like that of a leaky integrate-and-fire model, then, in contrast to the earlier reports that did not include the voltage effects, mutual excitation could lead to stable synchrony provided the synaptic reversal potential is below the voltage level the neuron would have reached in the absence of the interaction and threshold reset. This level is controlled by the applied current and the leakage parameters. Such synchrony is contingent on significant phase response (that would result, for example, by a sharp PRC jump) occurring during the synaptic rising phase. The rising phase, however, does not contribute significantly if it occurs before the voltage spike reaches its peak. Then a stable near-synchronous state can still exist between type 1 PRC neurons if the PRC shows a left skewness in its shape. These results are examined comprehensively using perfect integrate-and-fire, leaky integrate-and-fire, and skewed PRC shapes under the assumption of the weakly coupled oscillator theory applied to synaptically coupled neuron models.

1  Introduction

Oscillating neurons interacting with synaptic connections abound in a number of brain nuclei (Jahnsen & Llinás, 1984; Bevan & Wilson, 1999; Bolam, Hanley, Booth, & Bevan, 2000; Hutcheon & Yarom, 2000). Many of these connections are excitatory (Bevan, Magill, Terman, Bolam, & Wilson, 2002), and some are inhibitory (Filion & Tremblay, 1991). Some amount of coherent spiking (which can be termed “synchrony” or “transient synchrony”) (Wang & Rinzel, 1992; Bergman et al., 1998) (Salinas & Sejnowski, 2000; Engel, Fries, & Singer, 2001; Buzsáki, 2006) is a common phenomenon among such oscillating neurons. Tracing such coherent behavior to the mechanism of connectivity between the neurons and to the external synaptic driving influenced by the network (Stroeve & Gielen, 2001; Terman, Rubin, Yew, & Wilson, 2002; Lewis & Rinzel, 2003; Galán, Fourcaud-Trocmé, Ermentrout, & Urban, 2006; Tsodyks, Mitkov, & Sompolinsky, 1993) is crucial in order to systematically understand such behavior and control it for the purpose of therapeutic intervention. For example, networks in the basal ganglia in animal models of Parkinson's disease were observed to display spike time coherence (Raz, Vaadia, & Bergman, 2000; Wichmann, Bergman, & DeLong, 1994). Desynchronzing neuronal activity by deep brain stimulation (Wichmann & Delong, 2006; Wilson, 2013) is an effective therapeutic mechanism. These networks are complex, and often emulation of the behavior by equally complex network simulations (Guo, Rubin, McIntyre, Vitek, & Terman, 2008), while being helpful in engineering the required responses, still leaves unanswered several basic questions such as the relation between the intrinsic response properties of the neurons and the network behavior. For this, one must develop sufficient understanding of the relationship in pairs of coupled neurons and then extend it to include the network structure.

The intrinsic response behavior of an oscillating neuron is best captured by its phase response curve (PRC), which is easily measured experimentally (Winfree, 2001; Tateno & Robinson, 2007; Gouwens et al., 2010; Netoff, Schwemmer, & Lewis, 2012). (However, briefer and weaker stimuli will have be to employed in measuring such responses, or iPRCs, for the results of weakly coupled theory that are presented here to be applicable.) Thus, relating, the PRC to network behavior is insightful and rewarding in understanding neuronal synchrony, particularly when the neurons are self-oscillatory (Ermentrout, Galán, & Urban, 2007; Achuthan & Canavier, 2009; Smeal, Ermentrout, & White, 2010; Krogh-Madsen, Butera, Ermentrout, & Glass, 2012; Miranda-Domínguez & Netoff, 2013). However, real PRCs occur in a variety of shapes (Tateno & Robinson, 2007; Schultheiss, Edgerton, & Jaeger, 2010), although for ease of analytical studies, they may be broadly categorized into just two kinds, type 1 and type 2 (Hansel, Mato, & Meunier, 1995; Ermentrout, 1996), depending on whether the phase response involves for most of the neuron oscillation an advancement of the phase. However, some type 1 PRC neuronal models do show a negative, albeit small in amplitude, phase response immediately following the spike peak. In this work, we selected for analysis models that show purely positive PRC response, and a conductance-based neuronal model that has a type 1 PRC but with a brief negative phase response for simulating large coupled neurons. Considerable progress was made in relating type 1 and type 2 PRCs to network behavior (Van Vreeswijk, Abbott, & Ermentrout, 1994; Hansel et al., 1995; Ermentrout, 1996). The leaky integrate-and-fire (LIF) model (Gerstner & Kistler, 2002) and the perfect integrate-and-fire (PIF) are the simplest examples of type 1 PRC neuron model; the Hodgkin-Huxley model (Hodgkin & Huxley, 1952; Hansel, Mato, & Meunier, 1993) is an example of type 2 PRC neuron models. These models have been widely used as representative examples of these two classes of PRCs. However, the LIF and PIF models have sudden rises in the PRC level near zero phase; the LIF model also has an exponential PRC profile, and the PIF model results in a constant response throughout the oscillation. The HH model shows extremely little response near zero phase, a behavior quite common in a number of other models where certain currents can make the phase response shallow for a considerable portion of the oscillation following a spike. These features warrant a careful examination of the role of the PRC shape on the network behavior with the goal of discovering the effect of such features on the network synchrony.

Relating the PRC types (type 1 and type 2) to their network responses (whether synchrony ensues or fails) is not a simple task due to the recognition that the shape of each of these types can not only be modulated easily by the presence of certain currents, but can result in markedly different dynamical properties than their type would suggest when allowed to interact with other neurons in a network. Arthur, Burton, and Ermentrout (2013) showed that a theta model with a symmetric type 1 PRC, while retaining its oscillation period, can modulate its PRC shape when a small adaptation current is added to it. In particular, a skewness is effected in its shape. Such shape changes can significantly affect the neuron's interaction with other networked neurons. Delayed rectifier and the fast sodium currents can sensitively modulate the PRC shape and can cause PRC shape skewness (Ermentrout, Beverlin, & Netoff, 2012). In model computations, it was also shown that the position of the PRC peak can be sensitively affected by the presence of a persistent sodium current (Pfeuty, Mato, Golomb, & Hansel, 2005).

Here we selected three models—leaky integrate-and-fire (LIF), perfect integrate-and-fire (PIF), and a PRC model whose shape is altered using a parameter—to systematically study the relation between the PRC shape and the ensuing synchrony between a pair of coupled oscillating neurons using the approach of weakly coupled oscillator theory. The goal of this work is to not only exhaustively study these models such that prior results could be understood in the larger context of parameter spaces, but also to address how PRC shape skewness can alter the coherent behavior of coupled model neurons. Thus our results complement those reported earlier on LIF and related models, in particular those in Van Vreeswijk et al. (1994), Hansel et al. (1995), and Ermentrout (1996) and advance the role of PRC skewness on the synchrony mechanism in a variety of networks such as those reported in Pfeuty et al. (2005) and G. B. Ermentrout et al. (2012).

The LIF and the PIF neuron networks are some of the most popular networks being studied in neuroscience owing to their simplicity in formulation and their predictive capacity to other networks comprising more complex neuronal models (Lapicque, 1907; Van Vreeswijk et al., 1994; Hansel et al., 1995; Lewis & Rinzel, 2003; Jolivet, Lewis, & Gerstner, 2004; Burkitt, 2006; Ladenbauer, Augustin, Shiau, & Obermayer, 2012). Although they were not viewed from the point of view of the PRCs at the time of their introduction, they now serve as the simplest models that show drastic differences in their PRC shapes. The LIF model is more than a century old (Lapicque, 1907) and was first proposed by an experimenter to formulate accumulation of charge on the cell membrane in response to an applied current. Despite the dynamics of only the leakage current in the model, it nevertheless became a test bed for studies on synchrony and activity of neurons (Brunel & Sergi, 1998; Kuhlmann, Burkitt, Paolini, & Clark, 2002; Feng, 2003; Brette & Gerstner, 2005). With the realization of the importance of the PRC or phase response function in determining the synchrony of a network of oscillating neurons (Winfree, 1977; Tsubo, Teramae, & Fukai, 2007; Perez Velazquez et al., 2007; Tateno & Robinson, 2007; Achuthan & Canavier, 2009; Cui, Canavier Carmen, & Butera Robert, 2009; Smeal et al., 2010) the LIF model became more useful in constructing theories that may be extended or generalized to other more complex models (Chow & Kopell, 2000; Pfeuty, Mato, Golomb, & Hansel, 2003).

The work we present attempts a systematic study of the role of the PRC shape in ensuing stable synchrony or antisynchrony of mutually coupled neuron models. The work is mainly analytical with numerical results augmenting the analytical results. We first present the results for simple systems by illustrating them with numerical simulations and then introduce more complex phase response curves, along with simulation of large number of coupled neuron models. In section 2, the general phase reduction method of weakly coupled oscillator theory is briefly reviewed, along with the method of finding the stability of synchrony and antisynchrony. As a first simple example, the phase-independent PRC shape is studied in detail in section 3 by invoking the PRC of the PIF model. The jump near zero, together with a rising PRC profile, is studied in section  4 by examining the LIF model network of a pair of coupled neurons. This is again illustrated extensively with numerical simulations, and the analytical results delineating the synchrony and antisynchrony regimes are illustrated in several parameter spaces. In both models, the voltage time course is included in the analysis, and it plays a major role in the stability of synchrony, particularly when the coupling is excitatory. Also, the synaptic conductance used is always initiated at zero phase, thus having rising phases when the PRC is significantly nonzero. In section 5, we introduce more complex PRC shapes by first simulating a set of 100 mutually excitatory Wang-Buzsáki model neurons and then using PRC shapes constructed to deviate from the canonical PRC.

2  Ermentrout's Phase Reduction Method

Our goal is to determine the steady states of the coupled neurons and their linear stability. In addition to numerical simulations, we employ the weakly coupled oscillator theory that involves phase reduction to study the stability of synchrony and antisynchrony. The phase reduction method has been extensively used in neuroscience and applied mathematical studies, where in oscillators attracted to their limit-cycles receiving perturbation by either external stimulus or coupling can be reduced to a set of phase coupled oscillators (Winfree, 1967; Kuramoto, 1984; Ermentrout & Kopell, 1986; Brown, Moehlis, & Holmes, 2004). However, this is not the phase computed in Cartesian coordinates, but a nonlinear function representing the time shift the original oscillator would have to undergo to accommodate the perturbation. If the oscillator's variable is , then the perturbed variable is written as , where is the slow timescale, is the strength of perturbation, and is the phase of our interest. This was originally introduced by Neu (1979a, 1979b) and then generalized and formalized by Ermentrout (1981; Ermentrout & Kopell, 1991) using Fredholm's alternative, and then later restated as Malkin's theorem by Hoppensteadt and Izhikevich (1997). An excellent review of this method is found in Lewis and Skinner (2012).

The phase evolution can be expressed in the slow timescale after averaging over the fast timescale of the oscillator. When applied to neural oscillators, a neuron regularly spiking with a period and whose voltage variable is interacting with a second neuron via synaptic coupling receives a coupling perturbation proportional to where is the synaptic reversal potential, is the synaptic conductance delivered by the coupled neuron, and is the coupling strength. Similarly, the second neuron with a voltage variable receiving an identical coupling receives a current proportional to , where is the synaptic conductance received from the first neuron. In the reduced phase coupled equations, the phase derivatives are proportional to the convolutions of these currents with the phase response curve, termed the “interaction function.” Thus, the reduced phase evolution equations for two identical and mutually coupled neurons are written in terms of their phases and as
formula
2.1
formula
2.2
The time in the above equations is the slow timescale that we alluded to, but the symbol in the definition of the interaction function defined below is the oscillator's time used here as an integration variable:
formula
2.3
Note that the changes in frequency must be small for the theory to be valid, and in the normalized units, small changes in frequency are equivalent to small changes in the period. The “periodized” conductance is defined below. is the PRC. The phase is zero when the voltage time course reaches the peak of the action potential. It reaches the peak again at the end of the oscillation period, . The PRC obtained in this method that uses weakly coupled oscillator theory is sometimes termed “infinitesimal PRC” (iPRC) to indicate that it is obtained by using infinitesimally small stimulus strength such that the oscillation period is not significantly altered and the integration limits in the above integral are still meaningful in the coupled state of the neurons. The term
formula
2.4
which we call the “shunted PRC,” assumes importance in the study of the stability of synchrony (Hansel et al., 1995). This represents the modifications to the phase response introduced by the effect of the voltage and the synaptic reversal potential. The shunted PRC can become monotonically increasing or decreasing in some simple cases, leading to its direct relationship with stability.
If the synaptic conductance does not completely decay to zero before the start of another spike event, the remnant conductance is added to the newly triggered conductance level. Thus, in the steady state, a profile of the synaptic conductance is obtained that is termed “periodized” synaptic conductance and is computed as in Van Vreeswijk et al. (1994) and Ermentrout and Terman (2010). For an alpha function synaptic conductance where is the inverse of the synaptic time constant, it is given by
formula
2.5
has units of frequency. It represents the contribution of all the previous synaptic events to the interaction occurring at time and is used in the analysis in place of . The conductance due to the contribution of the previous spikes causes to have its peak () earlier than the peak time of the pure alpha function, particularly at high frequencies. For example, for a cell firing below 50 Hz with ms, the peak continues to be very close to 3 ms, but for 100 Hz, it occurs at 2.63 ms, and for 500 Hz, it occurs at 0.89 ms. Taking the difference of the two equations, 2.1 and 2.2, the evolution equation for the phase difference can be written as
formula
2.6
where is the growth function. Since can be expressed as a sum of its odd part and its even part , the right-hand side of the above equation can also be written as . This notation is used by some authors to study the stability (Van Vreeswijk et al., 1994; Ermentrout, 1996).
Synchrony () and antisynchrony () are both the steady states of this equation as it can be easily verified from equation 2.6 that becomes zero at these two steady states. Other phase-locked states could also exist, but they must be computed numerically. To become useful, the steady states must be stable. To test the stability of the synchrony, we must examine the eigenvalue . This is the slope of the growth function at . The parameter region with the property is the region of stable synchrony, and that with is the region of unstable synchrony. The curve in the parameter spaces obtained by setting is the critical curve across which stability changes occur. From the relation between and in equation 2.6 and the expression in equation 2.3, we can directly write as
formula
2.7
The antisynchrony is represented by the solution , and its stability is given the slope of the growth function at ; that is, the eigenvalue is given by
formula
2.8
In the steady state, both neurons maintain a constant phase difference, , and spike with an identical rate: . This common rate might be slightly different from their intrinsic oscillation frequency (). By adding equations 2.1 and 2.2 in the steady state (at ), we obtain the frequency in the steady state (Hansel et al., 1993) as
formula
2.9
The weakly coupled oscillator theory is strictly valid in the limit of the frequency being close to the uncoupled frequency. Having a nonzero is not a violation of the weakly coupled theory. The magnitude of , which, in our definition, can be controlled by , must be much smaller than .

3  Effect of Phase Independence of the PRC on Synchrony

The perfect integrate-and-fire (PIF) neuron model displays a phase-independent PRC. It has a voltage variable [] that linearly grows from a reset value () to a threshold level (). If the voltage in its ongoing evolution crosses the threshold level, it is assumed to have emitted a spike and thus is assigned the reset level. The spike itself is not represented in the model. The model is defined by a single equation for the voltage whose time derivative is a constant: such that is reset to when crosses . is the external constant applied current that is used to trigger oscillations in the model. The dynamics of two such coupled identical neurons are shown in Figure 1 by solving the equations for the two coupled neurons whose voltages are represented by and :
formula
3.1
formula
3.2
where the synaptic conductances triggered by the spike events in the first and the second neuron are, respectively, and ; is the synaptic reversal potential; and is the constant maximum synaptic coupling strength. We let the synaptic conductances assume a normalized alpha-function profile such that the synaptic conductance followed by a spike event at is given by , where is the synaptic rate constant. Since the applied current is identical to both oscillators, the neurons are identical oscillators, differing only by their initial states. The simulations are shown for mutual excitation. The two neurons that started with a voltage difference of 0.3 evolved into a synchronous state at (see Figure 1a). But when the current is reduced by half, the same initial conditions led to an antisynchronous state (see Figure 1b, top) that coexists with the synchronous state (see Figure 1b, bottom), which is obtained by choosing a different set of initial states.
Figure 1:

Numerical simulation of a network of two “perfect” integrate-and-fire neuron models. Coupling is excitatory, but results are qualitatively similar when the coupling is inhibitory. (a) The two voltage time courses at approaching a synchronous state. (b) The same initial conditions as in panel a leading to antisynchronous state, but more closely spaced initial conditions leading to a synchronous state when is lowered. Synaptic conductance is an alpha function with . , . .

Figure 1:

Numerical simulation of a network of two “perfect” integrate-and-fire neuron models. Coupling is excitatory, but results are qualitatively similar when the coupling is inhibitory. (a) The two voltage time courses at approaching a synchronous state. (b) The same initial conditions as in panel a leading to antisynchronous state, but more closely spaced initial conditions leading to a synchronous state when is lowered. Synaptic conductance is an alpha function with . , . .

These results may appear to be counterintuitive at first because it is generally assumed that neurons of this kind, which show type 1 PRCs (Hansel et al., 1995), should fail to synchronize for mutual excitation. The PIF model does display a type 1, PRC. Indeed, the voltage evolution of the uncoupled PIF model is simply a linearly growing function of time,
formula
3.3
and the phase response curve is simply the inverse of the time derivative of the voltage, which becomes
formula
3.4
Thus, the PRC is not only positive throughout but is a constant and independent of the phase. Consequently, the shunted PRC, equation 2.4, is a linear function of time. For , the shunted PRC is a linearly decreasing function of time but is still positive throughout. Hence, we expect the interaction function itself, equation 2.3, to be positive. This indeed is the case and can be seen more precisely by computing . For the alpha-function synaptic conductance, we obtain the interaction function as
formula
3.5
where , , . The time period itself is the time elapsed between spike reset and spike threshold and is obtained from equation 3.3 as . The interaction function is plotted in Figure 2a for which , , and , and the positive shift in the frequency at zero phase is . Thus, in the steady state of zero phase difference (synchrony), the oscillators are uniformly sped up from to . And at the synchronous state, the interaction function grows linearly, as can also be seen by Taylor-expanding near , , indicating that if the second oscillator leads the first by a small phase, then the first oscillator increases its frequency slightly proportional to the slope of the at . The second oscillator slows down because the slope of is negative at . This implies that the synchrony becomes stable, as is also indicated by the negative slope of the growth function (see Figure 2b, curve for ) at .
Figure 2:

Analytical prediction of synchrony and antisynchrony between two perfect integrate-and-fire neuron models. (a) The interaction function and the alpha-function synaptic conductance at an that makes the coupling excitatory. The PRC of the model is simply . (b) The growth function at three different . , the phase-locked states and their stability, are independent of , and thus panels b–d remain unaltered between excitatory and inhibitory couplings. (c) The steady phase differences as a function of the frequency displaying stable synchrony (filled circles resembling thick lines at ordinate values of 0 and ) at all frequencies, stable antisynchrony (filled circles at ) at low frequencies. The unfilled circles represent unstable steady-state phase differences. (d) Parameter display of synchrony (shaded, filling all the space) and antisynchrony (hatched) in the space of synaptic time constant () and firing rate (). Unit of ms is assigned to time for convenience.

Figure 2:

Analytical prediction of synchrony and antisynchrony between two perfect integrate-and-fire neuron models. (a) The interaction function and the alpha-function synaptic conductance at an that makes the coupling excitatory. The PRC of the model is simply . (b) The growth function at three different . , the phase-locked states and their stability, are independent of , and thus panels b–d remain unaltered between excitatory and inhibitory couplings. (c) The steady phase differences as a function of the frequency displaying stable synchrony (filled circles resembling thick lines at ordinate values of 0 and ) at all frequencies, stable antisynchrony (filled circles at ) at low frequencies. The unfilled circles represent unstable steady-state phase differences. (d) Parameter display of synchrony (shaded, filling all the space) and antisynchrony (hatched) in the space of synaptic time constant () and firing rate (). Unit of ms is assigned to time for convenience.

The expression for , equation 3.5, reveals that its dependence on is independent of , which merely alters the frequency of neurons uniformly. So for excitatory coupling, the frequency is elevated, and for inhibitory coupling, it is lowered, but the steady states and their stability are unaltered. This consequence is entirely due to the PRC being constant throughout the phase. The phase dependence of is entirely due to the voltage and . At low rates such that decays completely to zero before the onset of the next spike, becomes equivalent to the alpha function itself, and computed for slightly nonzero will have the contribution from the that is now mostly acting prior to the next spike time where is lowest. This leads to a rapidly decaying with , but as is advanced such that the influence of near the zero phase is minimized, the begins to increase linearly with . This leads to skewed , which will result in a nonzero phase-locked state. But at very high frequencies such that has no noticeable decaying tail, becomes more sinusoidal with dominating first Fourier mode and is similar to the well-known Kuramoto phase model. At intermediate firing rates, a stability switch can occur for the antisynchronous state when the synaptic conductance at half period is equal to inverse of the period. This can be easily seen when the stability of the antisynchrony is computed. The growth function defined in equation 2.6 is identical to
formula
3.6
and can be computed using equation 3.5 as
formula
3.7
which is independent of the reversal potential. And thus, as is already noted from , the type of synaptic coupling (whether excitatory or inhibitory) does not dictate the stability of the steady states. The growth function is plotted in Figure 2b as the rate is altered. At very low rates, the negative slope of at zero phase is not easily discernible, but its stability can be seen in the bifurcation diagram as a function of rate displayed in Figure 2c. In the next section, we show more systematically that synchrony is a stable solution at all rates by computing its eigenvalue and also show that the antisynchrony becomes a stable solution when . These results are summarized in Figure 2d.

3.1  Stability of Synchrony and Antisynchrony

3.1.1  Synchrony

Using the expressions in equations 3.3, 3.4, and 2.5 in equation 2.7, we get the eigenvalue for synchrony as
formula
3.8
This can be easily shown to be negative for .1 By noting that is a hyperbolic function that is purely positive and its Taylor series can be reexpressed as where for , we find that
formula
implying that the synchrony is stable. Using the method by Hansel et al. (1995), we can see how this negative sign arises due to the dominant negative contribution of the conductance downstroke overpowering the positive contribution due to the rising phase. For this, we first note that can be computed directly from the integral expression for in equation 2.3 via the relationship between the growth function and given in equation 2.6. By substituting only the expressions for and and not , we obtain
formula
3.9
In obtaining this, we used the fact that the integral of the derivative of the periodic function (i.e., ) is zero as a consequence of the fundamental theorem of calculus; this is also used below in computing the integrals. Since the slope of is positive before it reaches its peak at and is negative afterward, we can conveniently break up the above integral into two parts, each having a segment of that does not change the sign of its slope between the integration limits:
formula
3.10
Just as the mean of a function computed using a probability density function lies within the integration limits, the integral becomes equal to the product of and where , and is a fixed value between (and including) 0 and . This is guaranteed by the mean value theorem (see Courant, 2011). Similarly, the integral becomes equal to or simply , where is a fixed value between (and including) and . Thus, we can write equation 3.10 as
formula
3.11
The sign is asserted because is a monotonically increasing function, and thus . The value of can become zero only if itself were zero. And and become identical (and zero at the same time) only when itself were zero. We can also easily see from equation 3.8 by expanding the term in Taylor series that in the limit of either or going to zero, becomes 0. Thus, unless or becomes zero, the eigenvalue is negative (i.e., ). Thus, the synchrony is always stable for oscillating neurons () interacting with a nonzero () synaptic rate constant. It is also clear that the reason for the stability is the dominant contribution of the conductance decay phase () over the rising phase ().

3.1.2  Antisynchrony

The stability of antisynchrony is determined by the eigenvalue computed at (i.e., the sign of ). The region where the condition is satisfied indicates stable antisynchrony, marks unstable antisynchrony, and the curve in the parameter spaces where is a critical curve. can be computed directly from the expression in equation 3.7 as
formula
This can also be written in terms of the synaptic conductance at half period, , as
formula
3.12
where can be computed from equation 2.5 by inserting for . Thus, antisynchrony is stable in the region given by
formula
3.13
As increases (i.e., as the rate decreases), the left-hand side ( at half period) and the right-hand side () of the relation in equation 3.13 both decrease, but the conductance decreases faster than the rate due to its exponential nature. Thus, the above stability relation is easier to satisfy at low rates leading to stable antisynchrony: the stable antisynchrony region falls on the left of the critical curve in Figure 2d.

3.1.3  Synchrony and Antisynchrony for Double Exponential Synaptic Function

Changing the alpha-function to a double exponential synaptic function, , where , makes the periodized synaptic function
formula
3.14
where and . The analysis for the stability of the synchronous state remains identical to that carried out above for the case of alpha-function synaptic conductance. Hence, the synchrony remains stable at all rates and time constants. Even the condition for antisynchrony remains identical to that in equation 3.13. In particular, the eigenvalue for antisynchrony can be computed as before and is expressed as
formula
3.15
where and are the rising and decay rates of the synaptic conductance. This again can be simplified to the relation identical to that in equation 3.12, and thus the condition for the stability of antisynchrony is again given by
formula
3.16
where is obtained from equation 3.14 by inserting for . This condition leads to a dependence of the stability boundary on the time constants. This dependence is illustrated in Figure 3. As before the stable antisynchrony region falls on the left of the critical curves in Figure 3. If the rise time of the synaptic conductance is fast, the peak of the conductance occurs early, leading to a reduced synaptic conductance at half period. Consequently can be made equal to only for sufficiently small (i.e., high rate). Thus, the critical boundary of antisynchrony moves to smaller or higher rates, as seen Figure 3.
Figure 3:

Analytical prediction of antisynchrony between two perfect integrate-and-fire neuron models using a double exponential synaptic function in the plane of firing rate and synaptic decay time. The region to the left of each curve (the hatched region) marks the stability of antisynchrony for the corresponding synaptic rise time. Synchrony is stable in the entire parameter space.

Figure 3:

Analytical prediction of antisynchrony between two perfect integrate-and-fire neuron models using a double exponential synaptic function in the plane of firing rate and synaptic decay time. The region to the left of each curve (the hatched region) marks the stability of antisynchrony for the corresponding synaptic rise time. Synchrony is stable in the entire parameter space.

4  Effect of PRC Jump Near Zero Phase

We invoke the leaky integrate-and-fire (LIF) model (Lapicque, 1907; Gerstner & Kistler, 2002) to examine the effect of the PRC jump at the zero phase. The model's PRC, unlike that of the PIF, does depend on the phase and in fact increases monotonically. Conventionally, the LIF model is studied by setting the reset to 0 and the threshold level to 1 for convenience. However, setting them to more realistic values would enable realistic synaptic reversal potential and synaptic time constants. So although the analytical results are identical for both of these choices, we present numerical results under both of them; there are some noticeable differences in the parameter regimes where stability transitions occur. For choosing a realistic voltage shape, we use the Wang-Buzsáki model. Under either interpretation, the LIF model displays the same type of PRC, which is of type 1 and has sharp jumps at either extreme of the curve.

4.1  Model

The total current [] flowing through an LIF model membrane whose capacitance is and voltage is , is equal to the leakage current in the absence of any other input. By applying a steady external current (or even by adjusting ), the model can be made to oscillate. This is achieved by introducing a threshold-level crossing with reset to . The spike itself is not represented by the model and is ignored in the analysis. Such a neuron in the presence of the applied current evolves with the time course
formula
4.1
where , and . The model oscillates with a period
formula
4.2
and its PRC is simply proportional to , which can be computed as
formula
4.3
where . When two such neurons (whose membrane voltages are , ) are coupled synaptically with a reversal potential , they interact with their conductances (, ), which we model here with alpha functions as in the previous section. In particular the periodized synaptic conductance has the same form as in equation 2.5. Thus, when a neuron reaches its spike threshold, the corresponding synaptic conductance is generated, and the interaction between two such identical neurons takes place according to the following two coupled equations:
formula
4.4
formula
4.5
The parameter determines whether the coupling is excitatory or inhibitory. Under the weak coupling assumption, the coupling conductance is the smallness parameter that would not significantly alter the period of oscillation of each neuron.

4.2  Simulations

In order to simulate these equations are generated self-consistently using two first-order differential equations (Ermentrout, 2002), and in the analysis we use the form of . Simulations of the above equations are shown in Figures 4 and 5 for two sets of parameters. The first set of parameters emulates a Wang-Buzsáki model's voltage time course, and and are set to realistic values. The results are shown in Figure 4. Here we consider mutually excitatory coupling and thus set mV and ms. At 100 Hz, the two coupled neurons gradually tend toward synchrony of spike times. At 50 Hz, two different sets of initial conditions exhibit two different phase differences, suggesting a bistability between synchrony and antisynchrony. At very low frequency, a state that is very close to synchrony or a near-synchronous state is realized.
Figure 4:

Simulation of a mutually excitatory network of two leaky integrate-and-fire neurons whose voltages are modeled according to those of the Wang-Buzsáki model. (a) The voltage time courses (, ) evolving in time and approaching a synchronous oscillatory state at a firing frequency of 100 Hz. Spikes are not part of the model, and hence are not drawn. On the right, the effect of the coupling on is shown; it is excitatory. (b) Time courses as in panel a but at 50 Hz. Two sets of initial conditions lead to two different states: antisynchrony for widely separated initial conditions and synchrony for closely separated initial conditions. Right: The difference of successive spike times of the neurons normalized to the instantaneous period of one of the neurons (i.e the phase difference) is shown for the two sets of initial conditions. (c) As in panel b but at 10 Hz. Synchrony and antisynchrony are both unstable, but a near synchronous state that is very close to the synchronous state is stable. For all the simulations, the parameters are mS/cm, mV, F/cm, mV, mV. mV, ms, and  mS/cm.

Figure 4:

Simulation of a mutually excitatory network of two leaky integrate-and-fire neurons whose voltages are modeled according to those of the Wang-Buzsáki model. (a) The voltage time courses (, ) evolving in time and approaching a synchronous oscillatory state at a firing frequency of 100 Hz. Spikes are not part of the model, and hence are not drawn. On the right, the effect of the coupling on is shown; it is excitatory. (b) Time courses as in panel a but at 50 Hz. Two sets of initial conditions lead to two different states: antisynchrony for widely separated initial conditions and synchrony for closely separated initial conditions. Right: The difference of successive spike times of the neurons normalized to the instantaneous period of one of the neurons (i.e the phase difference) is shown for the two sets of initial conditions. (c) As in panel b but at 10 Hz. Synchrony and antisynchrony are both unstable, but a near synchronous state that is very close to the synchronous state is stable. For all the simulations, the parameters are mS/cm, mV, F/cm, mV, mV. mV, ms, and  mS/cm.

Figure 5:

Simulation of a pair of mutually excitatory LIF model neurons with the conventional parameters without units: , , , , and . (a) Voltage time courses of the two neuron models exhibiting antisynchrony, nonzero phase-locked state, and synchrony at three different levels of the applied current. (b) The phase difference of the consecutive spike times of the illustrations in panel a displaying the evolution toward steady states. The other parameters are , which is above the threshold, and , which determines the shape of the alpha-function synaptic conductance. The membrane effect is included in the computation.

Figure 5:

Simulation of a pair of mutually excitatory LIF model neurons with the conventional parameters without units: , , , , and . (a) Voltage time courses of the two neuron models exhibiting antisynchrony, nonzero phase-locked state, and synchrony at three different levels of the applied current. (b) The phase difference of the consecutive spike times of the illustrations in panel a displaying the evolution toward steady states. The other parameters are , which is above the threshold, and , which determines the shape of the alpha-function synaptic conductance. The membrane effect is included in the computation.

In the second set of simulations, we use the conventional set of parameters of the LIF whose voltage oscillates between 0 and 1. Units are not assigned to the time or the voltage variables. Simulation results at three different levels of for the same set of initial conditions are shown in Figure 5 for excitatory coupling with and . At high frequency, the phase difference clearly becomes zero, indicating synchrony; at low frequencies, the antisynchrony is stable; and at intermediate frequency, a nonzero phase-locked state is stable.

Similar results can also be obtained by using a double exponential synaptic function as long as the synaptic rise time is not instantaneous. As we will see in the analysis below, the rising phase of the synaptic conductance contributes to the stability of synchrony for excitatory coupling (i.e., for ), and the decaying phase contributes to the instability. (This is counter to their contributions in the PIF network of the previous section.) Depending on how an exponential synaptic function is implemented (whether the rising branch of appears as a curve with infinite slope or not), it may be possible to realize stable synchrony even with exponential synaptic function. If the rising phase occurs within a refractory period during which the PRC is zero, then instability of synchrony can result. Also note that there are two factors that help in stability of synchrony. For the double exponential, the peak of the periodized synaptic conductance () (see equation 3.14) occurs later than that of the , and thus the negative slope persisting beyond the synaptic rise time is in favor of the stable synchrony. At higher rates, the fraction of the rise time in comparison to that of the decay time also increases, helping to decrease the destabilizing factor.

4.3  Interaction and Growth Functions

The reduction of the coupled system (see equations 4.4 and 4.5) to phase coupled oscillators is carried out as outlined in section 2, and the interaction function and the growth functions are computed using the relations in equations 2.3 and 2.6. Using the expressions for the voltage and the PRC, equations 4.1 and 4.3, the shunted PRC can be written as
formula
4.6
where the factor is given by
formula
4.7
Hansel et al. (1995) assumed that the shunted PRC is monotonically increasing. Under this assumption, the eigenvalue for stability of synchrony (see equation 2.7) becomes positive, making synchrony unstable. Here, our goal is to fully explore the limitations of this assumption. depends on the period of oscillation. Since the period is a logarithmic function of the parameters (see equation 4.2), it is not immediately clear how varies as a function of the rate and, consequently, how depends on the parameters. However, it is clear that if , then the term becomes monotonically increasing, supporting the assumption of Hansel et al. But if , it is monotonically decreasing, leading to a stable synchrony, as we will show below. We also see that the eigenvalue for synchrony itself is proportional to ; thus, near the transition to synchrony, the neurons possess long transients because the eigenvalue goes through zero.
Substituting for (from equation 4.6) and (from equation 4.1) in equation 2.3, we obtain the interaction function as
formula
4.8
where , , , , , and , , and . The growth function is written from the relation in equation 2.6 as
formula
4.9
For the parameters of the LIF model whose voltage is based on the Wang-Buzsáki model (see Figure 6a–c), two sets of growth functions are illustrated. For mutual excitation, the growth function (see Figure 6d) shows a negative slope at zero phase (stable synchrony) and positive slope at half period (unstable antisynchrony) at 100 Hz (A/cm). Here, . The antisynchrony branch retains its state until the frequency is lowered below about 60 Hz, where it becomes stable (see Figure 6e). The stability of antisynchrony is also controlled by other conditions that do not directly depend on . The synchrony branch remains stable until the frequency is lowered such that the monotonicity of is affected. When , becomes positive and the synchrony becomes unstable. The antisynchrony branch also changes its stability because its stability, as we will see, also depends on .
Figure 6:

Analytical prediction of the dynamics of a network of two coupled leaky integrate-and-fire neuron models whose voltages are shaped according to those of the Wang-Buzsáki model. (a) The voltage time course of the each uncoupled LIF model (thick lines) compared to that of Wang-Buzsáki model (thin lines). (b) The PRC of the LIF model. (c) The firing rate as a function of the applied current . (d, e) Mutual excitation. (d) Growth functions at three levels of . A/cm for 100 Hz, 1.7825 A/cm for 50 Hz, and 0.53 A/cm for 25 Hz. (e) Stable (filled circles) and unstable (unfilled circles) steady-state phase differences with frequency. A range of firing rates (or ) can result in stable synchrony. (f, g) Mutual inhibition. (f) Growth functions are illustrated at three frequencies. (g) Steady-state stability of synchrony and other states with .

Figure 6:

Analytical prediction of the dynamics of a network of two coupled leaky integrate-and-fire neuron models whose voltages are shaped according to those of the Wang-Buzsáki model. (a) The voltage time course of the each uncoupled LIF model (thick lines) compared to that of Wang-Buzsáki model (thin lines). (b) The PRC of the LIF model. (c) The firing rate as a function of the applied current . (d, e) Mutual excitation. (d) Growth functions at three levels of . A/cm for 100 Hz, 1.7825 A/cm for 50 Hz, and 0.53 A/cm for 25 Hz. (e) Stable (filled circles) and unstable (unfilled circles) steady-state phase differences with frequency. A range of firing rates (or ) can result in stable synchrony. (f, g) Mutual inhibition. (f) Growth functions are illustrated at three frequencies. (g) Steady-state stability of synchrony and other states with .

For mutual inhibition (see Figure 6f,g), the growth function retains its negative slope (and hence its stability) at all frequencies because by decreasing and thus causing it to assume negative values, and in particular below , remains negative, and thus the shunted PRC is monotonically decreasing. The antisynchrony displays synchrony only at low frequencies, as is also the case when the effect of the voltage time course on stability was not included (Lewis & Rinzel, 2003).

For the conventional parameters of the LIF model whose voltage ranges between 0 and 1 (see Figures 7a–7c), growth functions and bifurcation of stable and unstable branches for mutually excitatory network are shown in Figures 7d to 7g. As it happened for the parameters of Figure 6, the growth function has a negative slope at large (giving a large firing rate). This is counter to what has been believed thus far for the LIF model networks based on the assumption that the voltage plays little role in the coupling mechanism. At large , the PRC becomes less steep, and thus the time dependence of the factor dominates the monotonically increasing PRC to cause the shunted PRC to decrease monotonically (). This effect would not exist if we ignored the voltage effect. As the is decreased, however, the steepness of the PRC begins to dominate, and when , synchrony becomes unstable. The stability of antisynchrony is determined by as well as other independent factors. A nonzero stable phase-locked state merges with unstable synchrony at that is barely sufficient to evoke spiking. We recover the stability diagram as a function of the synaptic rate (see Figure 7f) as reported by Van Vreeswijk et al. (1994) at small enough such that , and the diagram resembles that of a mutual inhibitory network when is sufficiently high such that the shunted PRC is monotonically decreasing (see Figure 7g).
Figure 7:

Analytical prediction of the dynamics of a pair of mutually excitatory LIF model neurons that use the conventional parameters as those set in Figure 5. (a) Voltage time course of each uncoupled model neuron. (b) Phase response curve of each neuron. (c) Firing rate of each uncoupled neuron with steady input current. If time were measured in ms, then the firing rate would be measured in kHz, and achieving firing rates of the order of 10s of Hz would become an extremely sensitive function of . (d) Growth functions at , displaying stable synchrony as is increased across . (e) Steady-state phase differences with . Filled circles indicate stability and unfilled instability. (f) Steady-state phase differences with at . (g) Same as in panel f but with . The results in panel f are similar to those of Figure 1 of Van Vreeswijk et al. (1994).

Figure 7:

Analytical prediction of the dynamics of a pair of mutually excitatory LIF model neurons that use the conventional parameters as those set in Figure 5. (a) Voltage time course of each uncoupled model neuron. (b) Phase response curve of each neuron. (c) Firing rate of each uncoupled neuron with steady input current. If time were measured in ms, then the firing rate would be measured in kHz, and achieving firing rates of the order of 10s of Hz would become an extremely sensitive function of . (d) Growth functions at , displaying stable synchrony as is increased across . (e) Steady-state phase differences with . Filled circles indicate stability and unfilled instability. (f) Steady-state phase differences with at . (g) Same as in panel f but with . The results in panel f are similar to those of Figure 1 of Van Vreeswijk et al. (1994).

4.4  Stability of Synchrony and Antisynchrony

Using the same method as in the previous section that was first used by Hansel et al. (1995) we show that the LIF pairs display stable synchrony if where is given by the relation in equation 4.7. Since the denominator of must be positive to elicit oscillations, synchrony can occur without any further conditions as long as . This condition encompasses regimes of mutual excitation and mutual inhibition.

4.4.1  Synchrony

The eigenvalue for the stability of synchrony is computed as
formula
We split this integral into two parts such that each part contains the portion of that does not change its sign:
formula
where is the peak time of the periodized synaptic conductance . As in the previous section, by using the mean value theorem, the first integral can be written as a product of and , where is between (and including) 0 and , and . The second integral can be written as the product of and , where is between (and including) and . Thus, we write the eigenvalue as
formula
Since (for nonzero ), the bracketed term is always positive. Thus, the sign of is determined by the sign of :
formula
The critical curve of synchrony is thus given by , or
formula
and the region of stable synchrony is given by . Since the denominator of must be positive in order to elicit spiking, this condition results in the following condition for synchrony:
formula
4.10
The right-hand side of the above inequality is nothing but the level the voltage time course of the LIF model would have reached if not reset to after crossing (see equation 4.1). This does not mean, however, that the excitation is acting like inhibition to cause stable synchrony, because in the coupled state, the firing rate of the neuron increases (see Figure 4a, right). The synaptic coupling would indeed begin to act like an inhibitory synapse if is less than when synchrony becomes stable without any further conditions. We can see this more clearly if we rewrite (from equation 4.7) as follows in terms of the oscillation period:
formula
Since is always positive, from this form, it is easy to see that becomes negative (i.e., synchrony becomes stable) if (see Figure 8a, right) without further restrictions on the parameters. This is also when the coupling acts like an inhibitory connection. However, if , the stability condition implies that the term in the brackets must be negative or, in simplified form,
formula
4.11
Figure 8:

Regions of stability of synchrony and antisynchrony for the LIF model with the conventional parameters. (a) Stability region in versus plane as the reversal potential is decreased from being mutual excitatory to being mutual inhibitory. The shaded region is the region of stable synchrony, and the hatched region is the region of stable antisynchrony. Multiple regions of bistability, as well as nonzero phase-locked states (white space) also exist. (b) Stability regions in and space at different levels of . (c) Stability regions in and space at different levels of the applied current. The model parameters are as those in Figure 7.

Figure 8:

Regions of stability of synchrony and antisynchrony for the LIF model with the conventional parameters. (a) Stability region in versus plane as the reversal potential is decreased from being mutual excitatory to being mutual inhibitory. The shaded region is the region of stable synchrony, and the hatched region is the region of stable antisynchrony. Multiple regions of bistability, as well as nonzero phase-locked states (white space) also exist. (b) Stability regions in and space at different levels of . (c) Stability regions in and space at different levels of the applied current. The model parameters are as those in Figure 7.

4.4.2  Antisynchrony

Simplifications as above are not feasible for the antisynchronous state. But the eigenvalue can be found directly by taking the derivative of at , and after simplification, we obtain
formula
where
formula
4.12
and . Since and the exponential term are positive, the sign of is determined by and , and thus the critical regions of stability are bounded by the two curves and . A portion of the stable synchronous region (that was defined by ) also becomes a region of stable antisynchrony if, in that portion, the condition is satisfied. In summary, the stable antisynchrony regions are defined by either or both of (1) , and (2) , . is independent of , but because does depend on , changing alters the boundary of antisynchrony corresponding to the curve . itself could be plotted as a function of its parameters to verify its sign. (Large numerical values due to could be alleviated by using the positive multiplication factor .)

These analytical results are summarized in the parameter space of (, ) at varying levels of reversal potential in Figure 8a. As we have already seen in Figure 7e, there can be three distinct regions displaying a nonzero phase-locked state (whose phase difference increases with rate toward half period), antisynchrony, and synchrony as (or the rate) is increased. This pattern appears at large enough reversal potentials as in Figure 8a (left). However, as is decreased, the boundary marking the stability switch of the antisynchrony remains unaltered ( is independent of ), but the point of transition indicated by the arrow in Figure 7e moves to smaller . If the new transition point occurs before the boundary of antisynchrony, then that regime is marked by and and exhibits bistability (see Figure 8a, middle). At low enough (as directed by equation 4.10), the entire range of exhibits stable synchrony, and the regime of below the antisynchrony boundary exhibits bistability. This bifurcation diagram is similar to that for mutually inhibitory LIF model neurons even in the absence of the voltage effect (Lewis & Rinzel, 2003).

In Figure 8b, the results are presented in the parameter space of (, ), where the distinction between slow synapse (small ) and fast synapse (large ) becomes clearer. For very slow synaptic time constants (see Figure 8b, left), synchrony and antisynchrony exist in complementary regions, and the transition from synchrony to antisynchrony occurs when is increased above . For faster synaptic time constants, this boundary can be between either a bistable region and a nonzero phase-locked state (low firing rates), or between synchrony and antisynchrony as before (high firing rates) (see Figure 8b, middle, right). Finally in Figure 8c, the transitions are shown at different rates in the parameter space of (, ), which makes it clear that at high enough rates, most levels of , whether excitatory or inhibitory, lead to stable synchrony, a result that emerges entirely due to the effect of the voltage on the synaptic coupling. And a stable antisynchrony exists for excitatory coupling at slow synaptic rates.

5  Effect of PRC Shape Skewness and the Emergence of Near-Synchrony

Thus far have we analyzed simpler models that have completely positive PRC shapes. In this section, we first illustrate simulation results using a more realistic conductance-based model that has a skewed PRC (along with a negative phase response immediately after spike peak) and then introduce a special class of phase response curve with a shape parameter. This allows for selectively tuning the skewness of the PRC. When the PRC was phase independent or had a sudden rise at zero phase, synchrony depended on the synaptic rising phase taking advantage of the finite and large-positive PRC level during that phase. Alternatively, if the PRC has a jump at zero phase and the rising phase of the synaptic conductance falls after the zero phase, synchrony may still become stable at appropriate firing rates. Generating the synaptic conductance (Wang & Rinzel, 1992) by solving a voltage-dependent differential equation such as
formula
5.1
may not allow a significant length of the synaptic rising phase to occur after the spike peak. This may prevent stabilizing the synchronous state. Panels in Figures 9a to 9c illustrate this using the Wang-Buzsáki model (see the appendix for the model). The model with a small applied current (A/cm) and a scale factor that scales the rates of sodium inactivation and the potassium activation uniformly () shows oscillations at 10 Hz and displays a type 1 PRC that looks like a bell-shaped curve except for a small, insignificant negative regime near zero phase. But more important, it displays a sudden rise in the PRC level within the first 2 milliseconds after the spike peak. But the synaptic conductance modeled with equation 5.1 ( mV, mV, , ) has most of its rise time before the spike peak during which the PRC is zero, and hence the opportunity to enhance the stabilizing component to the eigenvalue for synchrony of the coupled neurons is lost. The growth function (see Figure 9b) has a positive slope at zero spike time difference. Failure of synchrony is also seen in a network of 100 coupled neurons whose spike times are displayed in Figure 9c. Introducing an appropriate synaptic delay (see the dashed curve in Figure 9a) brings the rising phase of the conductance to the regime of finite PRC level, but the growth function computed at these parameter levels still has a positive slope at the synchronous state (see Figure 9b). Two factors would help stabilize synchrony in this situation: (1) increase the frequency of oscillation (as in the LIF model, Figure 7e) by increasing the applied current while maintaining the sudden PRC rise near early phases, and (2) adjust the parameters of the model in such a way that the sudden PRC rise occurs even earlier than 2 milliseconds. However, in the Wang-Buzsáki model, simply raising makes the jump in the PRC shallower.
Figure 9:

Near-synchrony in the Wang-Buzsáki excitatory network. (a–c) Failure of synchrony despite the sharp rise of the PRC near small phases, or delay in conductance. , A/cm. (a) , the synaptic conductance of the Wang-Buzsáki model illustrated with and without delay, and PRC profile at small phases. (b) The growth function of two mutually excitatory neurons without and with ( ms) delay. (c) Spike times of a network of 100 coupled neurons without synaptic delay failing to exhibit synchrony. mS/cm, mV. The initial conditions are chosen randomly on the limit cycle orbit of the uncoupled model. (d–f) Spike times of the Wang-Buzsáki network as the skewness of the PRC is altered while keeping the firing rate at 10 Hz (by adjusting ) and . Corresponding to these three simulations, the PRCs are shown in panel g and the growth functions for two coupled oscillators are shown in panel h. (i) The stable (filled circle) and unstable (unfilled circle) phase-locked solutions as a function of . .

Figure 9:

Near-synchrony in the Wang-Buzsáki excitatory network. (a–c) Failure of synchrony despite the sharp rise of the PRC near small phases, or delay in conductance. , A/cm. (a) , the synaptic conductance of the Wang-Buzsáki model illustrated with and without delay, and PRC profile at small phases. (b) The growth function of two mutually excitatory neurons without and with ( ms) delay. (c) Spike times of a network of 100 coupled neurons without synaptic delay failing to exhibit synchrony. mS/cm, mV. The initial conditions are chosen randomly on the limit cycle orbit of the uncoupled model. (d–f) Spike times of the Wang-Buzsáki network as the skewness of the PRC is altered while keeping the firing rate at 10 Hz (by adjusting ) and . Corresponding to these three simulations, the PRCs are shown in panel g and the growth functions for two coupled oscillators are shown in panel h. (i) The stable (filled circle) and unstable (unfilled circle) phase-locked solutions as a function of . .

A near-synchrony among the coupled neurons may be realized by causing the PRC to have skewness such that its response near early phases is subdued with or without having the PRC peak level move toward longer phases. We study this aspect in the rest of the section and in Figures 9d-i, 9.10, and 9.11. As is decreased, the sodium inactivation and the potassium activation both become slower, and a network of 100 coupled neurons displays increased amount of coherence with decreasing , as can be seen in the increased level of spike time alignment in Figures 9d to 9f (the rate is kept constant at 10 Hz by adjusting the ). As is decreased, the negative feedback currents become slower, making the spike broader and causing the hyperpolarization to last longer. Consequently, a subdued response near early phases results (see Figure 9g); and the sudden initial rise in the PRC diminishes. A second feature of the model is that its PRC becomes asymmetric. It is shifted toward longer phases, although the PRC peak position is not affected substantially by . The effect of these two causes is the emergence of a stable nonzero phase-locked state in the vicinity of the synchronous state. This appears as the intersection of the growth function with the horizontal axis near zero with a negative slope. This is a dramatic change in the character of the shape of the growth function considering that it mostly resembles a sinusoidal curve for a canonical PRC shape (Ermentrout, 1996; the dashed curves in Figures 9g and 9h). and are the same as for the case of . The nonzero phase-locked state moves toward the synchronous state with decreasing as can be seen in Figure 9i.

5.1  Formulation of Skewed PRC

The emergence of the nonzero phase-locked state can be sensitively dependent on the fine PRC features such as the overall symmetry of the PRC, the subdued response at the early phases, and the shift of the PRC peak. To study the nuances of this PRC structure, we construct a PRC as below that allows its control by a parameter
formula
5.2
where is the oscillation period and is a nonnegative real number that causes the left skewness (i.e., a long tail toward left of the peak). is a positive coefficient dependent on introduced to control the PRC amplitude. If we set to 1, the PRCs decrease in their amplitude considerably as increases because of the algebraic power. Although would appear as a coefficient multiplying the interaction function, and hence its choice does not alter the steady states or their stability, it may be set to a numerically approximate level of so that the peak amplitude of the PRCs would be commensurate with that when . The PRCs are illustrated in Figure 10a. The above formulation provides modulation of the PRC between phases 0 and . Since the PRC is either zero or subdued near small phases, the rising phase of the synaptic function does not contribute significantly toward synchrony stability. Thus, the interaction function for a mutually excitatory network can be computed without using the voltage effect using the following formula:
Figure 10:

Effect of PRC skewness on synchrony. (a) PRC shapes given by equation 5.2 as is increased. (b) The growth function for two mutually excitatory neurons getting affected by the period and, in particular, the antisynchrony becoming unstable as is increased even when the PRC has a small skewness. (c) Change of stability of synchrony, antisynchrony, and the nonzero phase-locked states with the rate for excitatory network when . (d) Same as in panel c for an inhibitory network. (e, f) Similar to panel c and d, but for large skewness (), depicting only the movement at low rates. In panel e, the nonzero phase-locked state is in close proximity to the unstable synchronous state (which is weakly repelling) at low rates, making it a near-synchronous state. (g) The growth functions for excitatory network revealing the nonzero phase-locked state surging toward the synchronous state as the rate is lowered at large skewness (). (h) Movement of steady states of the excitatory network as the skewness is increased at 10 Hz. (i) Steady states displayed in versus space for the excitatory network.

Figure 10:

Effect of PRC skewness on synchrony. (a) PRC shapes given by equation 5.2 as is increased. (b) The growth function for two mutually excitatory neurons getting affected by the period and, in particular, the antisynchrony becoming unstable as is increased even when the PRC has a small skewness. (c) Change of stability of synchrony, antisynchrony, and the nonzero phase-locked states with the rate for excitatory network when . (d) Same as in panel c for an inhibitory network. (e, f) Similar to panel c and d, but for large skewness (), depicting only the movement at low rates. In panel e, the nonzero phase-locked state is in close proximity to the unstable synchronous state (which is weakly repelling) at low rates, making it a near-synchronous state. (g) The growth functions for excitatory network revealing the nonzero phase-locked state surging toward the synchronous state as the rate is lowered at large skewness (). (h) Movement of steady states of the excitatory network as the skewness is increased at 10 Hz. (i) Steady states displayed in versus space for the excitatory network.

formula
where is the periodized synaptic function. We have set the coupling coefficient to unity. For mutual inhibition, the coefficient becomes negative, and the corresponding interaction is the same as the above but with a negative sign.

5.2  PRC with No Skewness

When , becomes symmetrical and is identical to the canonical type 1 PRC studied by Ermentrout (1996) and Ermentrout & Terman (2010). We set . When the PRC is symmetric, the frequency does not alter the stability of the steady states. However, the rate can alter the eigenvalue of stability of the synchrony. For example, if the synaptic conductance is a simple exponential , the periodized function still retains the same shape . Then the symmetric PRC [] results in the following growth function for mutual excitation,
formula
where . This cannot show any fixed points other than synchrony and antisynchrony. But the slope of at is , which shows that not only the state of synchrony is unstable (because the slope is positive) but its repulsion becomes weaker at low rates. (The slope is proportional to at large .) At the same time, the antisynchrony is also less tightly held at longer : the slope of at is exactly the negative of that at . These changes help the emergence of stable near-synchronous state when there is skewness present in the PRC.

5.3  PRC with Small Skewness

When , acquires a small skewness (see Figure 10a). PRC skewness of the kind that we have introduced can induce a number of higher-order Fourier components owing to the subdued response near zero phase. Even after including many Fourier modes, one may not be able to account for the fine feature of the PRC near the zero phase. However, such a close fit is not necessary to predict the steady states and their stability accurately. Ninety-nine percent of the Fourier power (see Dodla & Wilson, 2013b, for a discussion on the power) is contained in the following three modes:
formula
where . The shape of the periodized synaptic function depends on the frequency of oscillation, and thus the resultant interaction function also changes its functional form with the rate; the resultant Fourier approximation will have at least four harmonics of both cosine and sinusoidal terms in order to accurately predict the steady states, resulting in little insight into the relation between the Fourier modes in the PRC and the shape of the interaction function. Thus, the interaction function may be computed either analytically or numerically from the original shape of the PRC. It is easily computable for and results in the following growth function for the exponential (and excitatory) synaptic function:
formula
where , , and are constants dependent on and : , , , and . In addition to the synchronous and antisynchronous states, this expression for the can show, at appropriate rates, two other stable nonzero phase-locked states, and . The stability of the synchronous state is given by the slope of at and is given by the following simplified expression:
formula
Since the factor is positive, is positive. Since decays more slowly than , the term inside the last pair of parentheses is also positive. Hence, , resulting in an unstable synchronous state at all rates. However, at large , the slope can be seen to be proportional to , whereas for the symmetric PRC, it was proportional to