## 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).

## 3 Effect of Phase Independence of the PRC on Synchrony

### 3.1 Stability of Synchrony and Antisynchrony

#### 3.1.1 Synchrony

^{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 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 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: 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 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

#### 3.1.3 Synchrony and Antisynchrony for Double Exponential Synaptic Function

## 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

### 4.2 Simulations

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

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).

### 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

#### 4.4.2 Antisynchrony

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

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.