Synaptic and Intrinsic Determinants of the Phase Resetting Curve for Weak Coupling

doi:10.1007/s10827-010-0264-1

J Comput Neurosci. Author manuscript; available in PMC 2012 April 1.

Published in final edited form as:

J Comput Neurosci. 2011 April; 30(2): 373–390.

Published online 2010 August 11. doi: 10.1007/s10827-010-0264-1

PMCID: PMC3059351

NIHMSID: NIHMS264901

Synaptic and Intrinsic Determinants of the Phase Resetting Curve for Weak Coupling

Srisairam Achuthan,¹ Robert J. Butera,³ and Carmen C. Canavier^1,²

¹Neuroscience Center of Excellence New Orleans, LA 70112, USA

²Dept. of Ophthalmology, LSU HealthSciences Center, New Orleans, LA 70112, USA

³Laboratory for Neuroengineering, School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, GA, 30332, USA

Corresponding Author: Srisairam Achuthan Neuroscience Center 2020 Gravier St., Suite D New Orleans LA 70112 ; Email: sachut/at/lsuhsc.edu 504-599-0915 phone 504-568-5801 fax

The publisher's final edited version of this article is available at J Comput Neurosci

See other articles in PMC that cite the published article.

Abstract

A phase resetting curve (PRC) keeps track of the extent to which a perturbation at a given phase advances or delays the next spike, and can be used to predict phase locking in networks of oscillators. The PRC can be estimated by convolving the waveform of the perturbation with the infinitesimal PRC (iPRC) under the assumption of weak coupling. The iPRC is often defined with respect to an infinitesimal current as z_i(ϕ), where ϕ is phase, but can also be defined with respect to an infinitesimal conductance change as z_g(ϕ). In this paper, we first show that the two approaches are equivalent. Coupling waveforms corresponding to synapses with different time courses sample z_g(ϕ) in predictably different ways. We show that for oscillators with Type I excitability, an anomalous region in z_g(ϕ) with opposite sign to that seen otherwise is often observed during an action potential. If the duration of the synaptic perturbation is such that it effectively samples this region, PRCs with both advances and delays can be observed despite Type I excitability. We also show that changing the duration of a perturbation so that it preferentially samples regions of stable or unstable slopes in z_g(ϕ) can stabilize or destabilize synchrony in a network with the corresponding dynamics.

INTRODUCTION

Transiently synchronized assemblies of neurons are believed to underlie cognitive functions (Buzsaki, 2006). In humans, synchronous firing activity between distant brain regions is also thought to be the basis of complex phenomena related to memory such as word recollection (Fell et al., 2001) and facial recognition (Rodriguez et al., 1999). Abnormal synchronization has been demonstrated in cognitive and movement disorders, such as Parkinson's disease (Hutchison et al., 2004), epilepsy (Huguenard and McCormick, 2007; Stelt et al., 2004; Traub and Jefferys, 1994), tremor (Hammond et al., 2007) and schizophrenia (Timofeev and Steriade, 2004; Uhlhaas and Singer, 2006).

What are the conditions under which synchrony among coupled neural oscillators can be achieved? Phase resetting theory is often employed to address this question (Achuthan and Canavier, 2009; Acker et al., 2003; Brown et al., 2004; Canavier et al., 1997, 1999; Ermentrout, 1996; Goel and Ermentrout, 2002; Gutkin et al., 2005; Hansel et al., 1995; Netoff et al., 2005a, b; Oprisan et al., 2004; Sieling et al., 2009; Winfree, 2001) under the assumptions that the neural oscillators are limit cycle oscillators, and the perturbations affect only the phase and not the amplitude of the oscillation; this simplification results in phase oscillators that are completely characterized by a single variable, the phase. A phase response curve (PRC) plots the normalized change in cycle period as a function of the phase at which a perturbation is applied. The phase is defined modulo 1, such that it is zero at the beginning of a cycle and 1 at the end. Mathematicians interpret a shortening of the period as an increase in the phase that advances the end of the cycle, and a lengthening as a decrease in the phase that delays the end of the cycle.

Here we address weakly coupled phase oscillators, and make an argument for characterizing the infinitesimal response of noisy biological oscillators to conductance as well as to current. Weak coupling theory (Ermentrout 2002; Ermentrout and Kopell 1990, 1991) makes use of the PRC to an infinitesimal input z(ϕ). Since z(ϕ) is usually obtained from real neurons via deconvolution of x(ϕ) with the PRC measured in response to x(ϕ) applied at various phases within the cycle, the same experimental data can be analyzed in terms of either the infinitesimal response z_i(ϕ) to a current pulse or z_g(ϕ) to a conductance pulse. Previously, deconvolutions have been performed both ways. Netoff et al., (2005a) obtained z_i(ϕ) whereas Preyer and Butera (2005) obtained z_g(ϕ). In this paper, we first show that the two approaches are equivalent, but that in some cases z_g(ϕ) provides more useful insights about the expected shape of the phase resetting curve in response to realistic perturbations. We also show that for neural oscillators, Type I excitability (characterized by a gradual transition from quiescence to repetitive spiking as the applied depolarizing current is increased) bifurcation type does not always imply a monophasic (Type I) phase resetting curve (Ermentrout 1996), and that this result can also be explained in terms of how z_g(ϕ) is sampled. Thus one cannot infer the bifurcation mechanism by which spiking arises directly from the PRC shape, because the characteristic Type I shape is only preserved near the bifurcation, and even then only for sufficiently slow synaptic time courses. We next find that the synaptic input conductance waveforms with a brief time course sample z_g(ϕ) more effectively in comparison to those with a longer time course. This explains the role of synaptic rise time and duration in determining synchronization properties (van Vreeswijk et al., 1994) of a network of two weakly coupled oscillators.

METHODS

Phase Model

The original formulation of the phase model due to Winfree (1967) defined z(ϕ) as a “sensitivity function”:

(1)

where Δω is the change in frequency due a perturbation S. The phase (ϕ) advances uniformly in time except for a perturbation in frequency equal to x(ϕ)z(ϕ) where x(ϕ) is the “influence function”. This can be written as follows:

(2)

with ε very small, representing weak coupling between oscillators. In general, ε will be considered to be incorporated within x(ϕ) except when we utilize ε explicitly in the section Weak Coupling Approach to Study Synchronization among Coupled Oscillators. The subscript j instead of i is meant to indicate that in a coupled system, the influence function depends upon the phase of the other oscillator(s), but for the purpose of understanding the phase resetting in an isolated oscillator, x(ϕ) can be parameterized in terms of the phase in the oscillator under study. Note that Winfree defined the sensitivity function in terms of the effect on the frequency, whereas we have referred to z(ϕ) as the infinitesimal phase resetting curve, which measures changes in normalized period (that is, phase) and not frequency. For weak coupling, the phase resetting (Δϕ) is quite small and can be considered equivalent to the normalized change in frequency (Δω). If the unperturbed period is T, then the perturbed period T_p is T(1−Δϕ), where Δϕ = (T−T_p)/T. It follows that the unperturbed frequency is ω=1/T, and the perturbed frequency = 1/(T(1−Δϕ)), so the change in frequency is actually Δω = 1/(T(1−Δϕ)) − 1/T. However, for normalized frequency T=ω=1 and small resetting as would be expected for weak coupling, this approach is approximately correct because Δω = Δϕ / (1−Δϕ) [congruent with]

Δϕ as Δϕ goes to zero (as in eqn. 3.7 in Hansel et al., 1995). In this study, we will follow the weak coupling assumption that the sensitivity function z(ϕ) is defined in terms of the phase resetting, or change in period, rather than in terms of frequency as in Winfree's original definition.

Direct Method to Measure Infinitesimal Phase Resetting Curve

There are two ways to measure the infinitesimal phase resetting curve: one is direct, and one involves deconvolution. The direct method uses a square pulse rather than a delta function because a delta function perturbation in current has no duration and theoretically produces an instantaneous change in membrane potential; such a discrete, discontinuous jump cannot be produced in a real-time, continuous system. We must obtain the response to a pulse and consider the limit as the area (not the amplitude) goes to zero (not one) in order to guarantee that we are measuring the response to a weak input so that we can assume linearity. In practice, the width of the pulse Δt must go to zero, otherwise the pulse is not delivered precisely at a phase of ϕ, but is spread out and some portion of the input is applied at later phases. The test pulse p(t) is applied at ϕ = t/T.

(3)

with a representing the amplitude of the pulse and u(t) the unit step function. The infinitesimal PRC z(ϕ) is then defined as follows:

(4)

Both the perturbation and the evoked phase resetting go to zero in the limit, but their ratio converges to the theoretical phase resetting in response to a delta function. Convergence is achieved in a scaling region in which halving the perturbation halves the response, but in practice when the perturbation becomes too small this region is lost to noise, even for models due to round off errors. If the perturbation is defined in terms of an instantaneous change in membrane potential (Brown et al., 2004), then it is accurate to ignore the duration of the pulse and let the amplitude go to zero (in fact this is how the adjoint is defined, see below) but this is not a perturbation that can be implemented physically, only mathematically. Frequently a current pulse (Galan et al., 2005; Gutkin et al., 2005) is applied to probe phase resetting, but a conductance pulse (Maran and Canavier 2008) can also be applied. If the phase resetting produced by a square pulse is normalized by the unit area of the pulse, the scaled quantity z(ϕ) becomes formally equivalent to the instantaneous phase resetting produced by a unit impulse, or delta function applied at phase ϕ. The subscript i or g indicates whether the amplitude of the square test pulses were measured in terms of current (i) or conductance (g). The accuracy of the measurement of the phase resetting places a lower limit on how far the amplitude of the pulse can be reduced prior to convergence to the true value of z(ϕ). Note that if the usual neurophysiological conventions for current are followed, a minus sign is required to convert from the sign convention for an applied current to a membrane current.

(5)

z_i(ϕ) by convention refers to the phase resetting induced by a depolarizing pulse of unit area. The response to a hyperpolarizing pulse is the same except the sign is inverted. Since iΔt has units of charge, z_i(ϕ) has units of phase change per unit charge. In the figures in this manuscript, we use charge densities instead of charge since the model parameters (see Appendices (A) and (B)) are given per unit membrane area.

Convolution and Deconvolution Methods

As noted in the Introduction, experimental studies have attempted to measure the infinitesimal phase resetting indirectly, via deconvolution; Netoff et al., (2005a) obtained z_i(ϕ) whereas Preyer and Butera (2005) obtained z_g(ϕ). This indirect method is possible if we assume that the phase resetting H(ϕ) due to an arbitrary input can be approximated by convolving the infinitesimal phase response z(ϕ) with the arbitrary input x(ϕ) as illustrated in Appendix (C) (see Fig. 6B); conversely, z(ϕ) can be approximated by deconvolving H(ϕ) and x(ϕ), although the presence of noise in real systems renders this a difficult problem (Khoo, 2000; Netoff et al., 2005a). The influence x(ϕ) and sensitivity z(ϕ) functions can be defined in more than one way. A common choice is to use z_i(ϕ) as defined above as the infinitesimal phase resetting curve with respect to current. The mathematical advantage of using z_i(t) is that it can also be calculated as dϕ/dV from the differential equations that describe the oscillator using the adjoint linear system (Brown et al., 2004; Ermentrout, 2002; Ermentrout and Kopell, 1991; Hoppensteadt and Izhikevich, 1997). The adjoint z_v(ϕ) is obtained in units of inverse membrane potential and should be scaled by the membrane capacitance C_m, in order to obtain z_i(ϕ) since iΔt = C_m ΔV. z_i(ϕ) also has the advantage of being quite general in that it reflects only intrinsic dynamics and is independent of the nature of the perturbation (i.e., excitatory or inhibitory). We can define the phase resetting H(ϕ) in terms of z_i(ϕ) as the circular convolution integral as follows:

(6)

where I_syn(ϕ) represents the synaptic current waveform, g_syn(ϕ) represents the synaptic conductance waveform, V_post(ϕ) is the postsynaptic membrane potential and V_syn is the synaptic reversal potential. The term (V_post(ϕ)−V_syn) is problematic because it causes the synaptic current to be a function not only of the conductance in the presynaptic (driving) oscillator but also of the voltage in the postsynaptic (driven) oscillator. Logically, we have included (V_post(ϕ)−V_syn) in the influence function (see Appendix of Oprisan and Canavier, 2002) because we treat the current as the external perturbation to the system. However, others (see eqn. 7.31 in Rinzel and Ermentrout, 1998; Netoff et al., 2005a; Ermentrout and Kopell, 1991) have grouped this term in the sensitivity function with z_i(ϕ) instead because both of these terms depend only upon the postsynaptic oscillator, whereas g_syn(ϕ) depends only on the presynaptic oscillator. Convolution is associative so the grouping does not affect the result. A priori knowledge about the voltage waveform of the driven system is required to scale the impulse response BEFORE it is convolved with the conductance waveform. In an experimental system, this waveform may be contaminated by noise and non-stationarity, thus it makes more sense to directly measure the phase response to a perturbation in conductance. This is made possible by the use of the dynamic clamp (Sharp et al., 1993a, 1993b), which computes the (V_post(ϕ)−V_syn) at each time point based on the direct measurement of the postsynaptic membrane potential.

FIGURE 6

Analogy between linear systems theory and phase resetting theory under the assumption of weak coupling. (A) The output y(t) of a linear time invariant system (LTI) produced by an input waveform x(t) is the convolution (represented as *) between x(t) and (more ...)

From an engineering point of view, the analogy with linear systems theory (see Appendix (C)) is improved if all aspects of the postsynaptic (driven) oscillator are included in the measurement of the infinitesimal phase response, and only the input external to the driven oscillator, which we define as the synaptic conductance waveform, is included in the arbitrary input. If we define the perturbation x(ϕ) with respect to the synaptic conductance waveform as in Preyer and Butera (2005) instead of the synaptic current, we have:

(7)

In the above equation z_g(ϕ) must be defined separately for each choice of synaptic reversal potential, hence it represents not only the intrinsic dynamics of the oscillator but also the interaction of the synaptic dynamics with the intrinsic dynamics. z_g(ϕ) is defined as follows:

(8)

where g is the amplitude of the conductance perturbation which is always non-negative. A dimensional analysis reveals that conductance can be measured in Farads per second. This is evident from i = gV, since current can be measured in Coulombs per second and potential can be measured as Coulombs per Farad as charge divided by capacitance gives the potential. Therefore, z_g(ϕ) has units of phase change per unit capacitance. In this case, it is not necessary to know (V_post(ϕ)−V_syn) a priori, therefore the arbitrary input consists simply of the conductance waveform that depends only on the external drive received by the oscillator under consideration.

In order to compute H(ϕ)= x(ϕ)*z(ϕ) (where the * denotes the convolution operation) directly from the infinitesimal PRC sampled at a finite number of phases, we need to use a convolution sum for discrete systems. In order to use continuous systems theory, one would need to assume an analytically tractable form for the PRC for which the integral is known, but this is an unnecessary and restrictive assumption when we can utilize the experimental data directly instead. x(ϕ) and z(ϕ) were sampled at 1000 evenly spaced values (n) throughout the cycle. The height of each pulse in the train was an estimate of the area under the curve (chapter 10 of Izhikevich, 2007) in the 1000 sampled increments. This was obtained by multiplying the sampled value of x(ϕ) by the interval (Δt) between pulses, which represents the total charge to convolve with the phase resetting per unit charge z_i(ϕ) or the total capacitance to convolve with the phase resetting per unit capacitance z_g(ϕ)as shown in eqn. (9) below.

(9)

To obtain H_i(ϕ), we replace in eqn. (9) z(ϕ) with z_i(ϕ) and x(ϕ−σ) with I_syn(ϕ−σ). On the other hand, to obtain H_g(ϕ), we replace in eqn. (9) z(ϕ) with z_g(ϕ) and x(ϕ−σ) with g_syn(ϕ−σ).

Weak Coupling Approach to Study Synchronization among Coupled Oscillators

Let ϕ₁(t) and ϕ₂(t) represent the absolute phases of two weakly coupled identical neural oscillators as a function of real time t. Note that these can be determined by integrating eqn. (2). Then the relative phase difference between them can be written as θ(t) = (ϕ₁(t)−ϕ₂(t)). From (2) we have:

(10)

Since ε is small due to the assumption of weak coupling between the oscillators, the rate of change of θ(t), i.e. dθ/dt, is slow in real time t. We can therefore average dθ/dt over a cycle period, T. Using a change of coordinates (Theorem 4.1.1 from Guckenheimer and Holmes, 1983),

(11)

gives us the critical simplification that the relative phase equation M12

can be considered to be approximately constant during the course of a single cycle so that we can perform the convolution integrals as follows:

(12)

For two identical phase locked oscillators, if one oscillator has a phase relative to the other of ϕ, the other has a phase modulo 1 relative to the first of −ϕ,. This implies that:

(13)

where

Note that the H function for each coupled oscillator is the same expression we obtained earlier via a different method to represent the phase resetting curve of each isolated oscillator. In eqns. (6) and (7), H(ϕ) represented the resetting of an oscillator due to an input perturbation applied at each phase, ϕ. Eqn. (13) involves the expression for H( equation M16

) which represents the phase resetting for the case of two oscillators coupled with relative phase, equation M17

. Also note that the H function appears in the derivative of the slowly varying relative phase (eqn. (13)) but not in that of the difference of the rapidly varying absolute phases (eqn. (10)) due to the contribution of averaging theory. Here the perturbation x(ϕ) is explicitly driven by the partner oscillator, and the two oscillators can lock at a common frequency when H(ϕ)=H(−ϕ) such that H(ϕ)−H(−ϕ) = 0. Let G(ϕ) = H(ϕ)−H (−ϕ). Linearizing the right hand side about a fixed point at a relative phase of ϕ using the first oscillator as the reference, the fixed point is a stable if G'(ϕ)<0 because a negative slope produces a decaying exponential (Ermentrout and Kopell, 1984; Van Vreeswijk et al.,1994). Therefore, the sign of slope of the odd portion of H(ϕ) determines the stability for weak coupling. Using the convention that phase advances produce a positive Δϕ, negative slopes are stable whereas positive slopes are unstable.

Simulations

The Morris-Lecar (1981) and the Wang and Buzsaki models are (1996) used as examples of simple neural models. The differential equations were integrated using a variable step size implicit fifth order Runge Kutta method (Hairer and Wanner, 1991). The model equations and parameters used in this paper are given in Appendices (A) and (B).

Direct Measurement of the PRC

PRCs were measured as in Acker et al. (2003). The presynaptic model neuron is initialized at its spiking threshold when the coupling is turned on for a single cycle of the presynaptic neuron at different phases in the cycle of the postsynaptic model neuron. Zero phase is defined as an upward crossing of the threshold potential defined below. The PRCs were generated for each neuron by stimulating a single repetitively firing neuron by driving the synaptic conductance in that neuron with a single action potential from an identical neuron. This type of PRC is often called a spike response curve (Acker et al., 2003). The single action potential in the presynaptic model neuron triggers a rise and fall in the postsynaptic conductance and thereby in the postsynaptic synaptic current, which serves as the perturbation x(ϕ) utilized to generate the PRC. The action potential threshold was set to a voltage at which the synaptic conductance trace began to noticeably increase above zero, which was −14 mV for a synaptic half activation (V_half) potential of 0.0 mV and 14 mV for V_half of 28 mV. We defined the start of the limit cycle (a phase of 0) to coincide with the time at which an input would be received from a presynaptic neuron firing in exact synchrony with the neuron receiving an input. This time occurs before the peak of the action potential due to the form of the coupling equations chosen. The first order phase resetting is defined as Δϕ₁ = (T−T₁)/T, where T₁ is the duration of the perturbed cycle containing the start of the perturbation. The second order resetting is the effect, if any, on the next cycle Δϕ₂ = (T−T₂)/T, and is often neglected.

RESULTS

Two Methods of Estimating the PRC to Realistic Perturbations are Equivalent

Figure 1A compares the iPRCs measured with respect to conductance (z_g(ϕ)), continuous gray curves) with those measured with respect to current (z_i(ϕ)(V_post(ϕ) −V_syn), dashed black curves) and scaled with respect to the driving force (V_post(ϕ) −V_syn). A Morris Lecar (ML) oscillator receiving an inhibitory synaptic input was used for this example. Eqns. (6) and (7) imply that:

(14)

for 0 ≤ σ ≤ 1, where 1 is the normalized period of the oscillating neuron. The essential difference between the two methods is that when the z_i(ϕ) is scaled, V_post(ϕ) is taken as the unperturbed waveform of the free running oscillator whereas when the z_g(ϕ) is measured, V_post(ϕ) takes on the actual values in the driven oscillator during the perturbation. In the limit of infinitesimally small perturbations, Fig. 1A shows that the two methods produce the same result.

FIGURE 1

Two ways of measuring the iPRC are equivalent. (A) z_g (gray curve) is compared to z_i scaled by (V_post(ϕ)−V_syn) (dashed curve) withV_syn = −75.0 mV corresponding to an inhibitory synapse for this example. (B) The directly measured (more ...)

For weak coupling, the PRC generated by a spike in the presynaptic neuron can be estimated¹ (Maran and Canavier, 2008, see also Netoff et al., 2005a; Rinzel and Ermentrout 1998; Preyer and Butera 2005) by convolving the perturbation in conductance or current with the appropriate infinitesimal PRC. Eqn. (13) implies that;

(15)

Figure 1B confirms that the PRCs generated by either of the convolution approaches for weak coupling approximate quite well the actual PRC measured by the direct method in the example given for the Morris Lecar model with weak resetting (only a small fraction of the cycle). See Appendix (D) for details on how the approximation can be improved.

Whether the PRC Contains both Advances and Delays Depends on z_g(ϕ) and Relative Coupling Duration

Frequently, the shape of the PRC is used to infer excitability type. Neurons exhibiting Type I excitability have a gradual onset of spiking with arbitrarily low frequencies as depolarizing current is applied from a quiescent state, whereas neurons exhibiting Type II excitability have an abrupt onset of spiking at a threshold level of current below which repetitive firing cannot be supported, (Hodgkin, 1948). PRCs that have only delays or advances are often categorized as Type I whereas those with a region of advances and a region of delays are categorized as Type II (Hansel et al., 1995). Since Type I excitability² is often associated with a saddle node bifurcation structure that makes it a function as an integrator, excitations usually produce advances and inhibitions usually produce delays, so Type I excitability usually implies Type I PRCs (Ermentrout, 1996). Similarly, Type II excitability is associated with a Hopf bifurcation and functions as a resonator. In this case, excitations produce early delays and late advances, and the converse applies for inhibition. These rules of thumb only strictly apply near the bifurcations. Proximity to a bifurcation can be defined in two ways. One is in terms of the trajectory in the state space: spiking neurons are generally close to the ghost of the saddle node during the interspike interval but not during the large voltage excursion experienced during the action potential. Thus Figs. 2A1 and 2B1 show that there is an anomalous region of advance in response to inhibition in z_g(ϕ) for a Morris Lecar Type I example at two frequencies. We use z_g(ϕ) instead of z_i(ϕ) because the multiplication by the driving force could in some cases alter the position of the anomalous region by flipping the sign. The presence of the anomalous region presents the possibility that despite the emergence of spiking via a saddle node bifurcation and Type I excitability, if the perturbation is short enough to preferentially sample this region, a Type II phase resetting curve could result. Figures 2A3 and 2B3 show that this is indeed the case. Furthermore, increasing the frequency preferentially shortens the interspike interval, and increases the fraction of the cycle occupied by the action potential, which in turn augments the anomalous region. Therefore, the same synaptic perturbation will sample the anomalous region more effectively at higher frequencies. Synaptic dynamics were held constant for the two frequencies, and the anomalous region is consequently much more prominent for the case with larger I_stim (15 μA/cm² and 13 ms period versus 9μA/cm² and a 26 ms period). The anomalous region in panels 2A3 and 2B3 is displaced to the left from that shown in panels 2A1 and 2B1 because of how the synaptic perturbation samples the anomalous region in z_g(ϕ). Purely Type I PRCs are only observed when the anomalous sign region is brief compared to the perturbation, so the region of anomalous sign is lost in the smoothing produced by the convolution.

FIGURE 2

For a constant synaptic input duration, higher frequencies resulting in an action potential that occupies a large fraction of cycle can lead to a Type II PRC with Type I excitability. (A) Low frequency example (I_stim = 9.0 μA/cm²). (A1) z_g(ϕ) (more ...)

The other way to define proximity is in terms of the bifurcation parameter, which is the depolarizing applied stimulus current I_stim. As I_stim is increased beyond the saddle node bifurcation at which two fixed points collide, spiking eventually terminates in either a sub (Fig. 2C) or supercritical Hopf. For the ML example, the SNIC bifurcation (Fig. 2C) occurs at I_stim = 8.326 μA/cm², and spiking is lost to depolarization block at the subcritical Hopf bifurcation at I_stim = 22.78 μA/cm². However, the dynamics at I_stim = 9 and 15 μA/cm² are still those of an integrator influenced by the SNIC, and not those of a resonator, thus the presence of an anomalous region in the PRC can persist far from the SNIC and does not appear to be due to the proximity to the Hopf. The Morris Lecar model has wide action potentials, and only two dimensions. In order to show the generality of this mechanism, we examine a second example in Wang and Buzsaki model neuron with thin action potentials, also for inhibition with Type I excitability, and also at two frequencies. The same principles apply, as z_g(ϕ) has a larger anomalous region for the higher frequency (Fig. 3A1 compared to Fig 3B1), resulting in a PRC with a more pronounced Type II character (Fig. 3A3 vs Fig. 3B3), although the anomalous region is less prominent than for the ML example in Fig. 2. For this example, the SNIC bifurcation (Fig. 3C) occurs at I_stim = 0.1601 μA/cm², and spiking is lost to depolarization block 17 at the supercritical Hopf bifurcation at I_stim = 25.13 μA/cm². However, the dynamics of the examples shown at I_stim = 2 and 8 μA/cm² are still those of an integrator influenced by the SNIC, thus again the results apply far from the bifurcation in parameter space.

FIGURE 3

Even models with thin action potentials can have Type II PRCs despite Type I excitability. (A) Low frequency example (I_stim = 2.0 μA/cm²). (A1) z_g(ϕ) has an anomalous region (above the dotted horizontal line). (A2) Membrane potential for (more ...)

The Slopes of the PRC Ddetermine the Stability of Phase Locking in Weak Coupling

We noted in the Methods that the slope of the odd portion G(ϕ) of the H function gives the stability of a fixed point. H(0) = H(1) for weak coupling since the H function is by definition periodic with period T or normalized period 1, and evaluating the slope of G(ϕ) at zero phase gives G'(0) = (G(Δϕ)−G(0))/Δϕ. Since G(0) = 0, G'(0) = (H(Δϕ)− H(1−Δϕ))/Δϕ = (ΔϕH('0) − (−ΔϕH'(0))/Δϕ = 2H'(0). This result generalizes to any fixed point because the common network period guarantees that H(ϕ) = H(1−ϕ) at the fixed point, which implies that the slope of the G function is the sum H'(ϕ) + H'(1−ϕ). When the coupling becomes strong, the cross product term H'(ϕ)H'(1−ϕ) must also be considered (Dror et al., 1999). For the purposes of this study, we focus on the result that the sign of the slope of phase resetting curve H'(ϕ) at zero is sufficient to give the stability of synchrony for identically coupled identical oscillators, which will be quite useful in developing an intuition for when synchrony can occur as we show in the next few sections.

The Slope of z_g Combined with Relative Coupling Duration Determines the Slope of the PRC

The slope of z_g(ϕ) but not z_i(ϕ) has important implications for the stability of phase locking, because multiplication by the driving force can change the sign of the slope. In addition to anomalous regions in the sign of z_g(ϕ) regions of anomalous slope can exist, such as the initial region of anomalous slope in the ML example in Fig. 4B. Whether this region is preserved in the PRC depends upon the level of filtering, or smoothing produced by the convolution, with longer duration synaptic interactions producing more smoothing than shorter ones. In this example, the frequency was kept constant, but two values for the time constant of synaptic decay were used, one fast (1 ms) and one slow (3 ms), to produce a brief perturbation and one with a longer duration. A brief perturbation (Fig. 4A1) samples z_g(ϕ) (Fig. 4B) in such a way that the PRC in Fig. 4A2 preserves the initial positive slope of z_g(ϕ) (inset to Fig. 4A2). The positive slope at a phase of zero implies that synchrony is unstable; therefore two ML neurons coupled via such synapses do not synchronize but instead fire in antiphase (Fig. 4A3). On the other hand, the longer perturbation (Fig 4E1) smooths the z_g(ϕ) by sampling a larger segment of z_g(ϕ), which causes the small anomalous initial region of positive slope in Fig. 4B to disappear in the corresponding PRC plot (Fig. 4E2). The smoothing effect is a result of the convolution in eqn. (7). As a consequence, two ML model neurons coupled with a longer duration can now synchronize (Fig. 4E3) since the slope of the PRC at phase zero is negative in Fig. 4E2, in contrast to the absence of synchrony with shorter duration coupling (Fig. 4A3). Thus, an examination of z_g(ϕ) explains how changing the duration of an inhibitory/excitatory synaptic input can affect the ability of neural oscillators to synchronize. In the weak coupling regime, the magnitude of the synaptic conductance strength does not change the shape of the PRC, hence the choice of conductance strength is arbitrary. We chose a smaller conductance (g_syn = 0.01 mS/cm²) for the slower time constant compared to that for the faster time constant (g_syn = 0.05 mS/cm²) in order to more closely balance the amount of resetting observed at each phase, although this balancing cannot be done systematically, because the ratio of the resetting produced by the short and long duration pulses itself is a function of phase. This also applies to Fig. 5.

FIGURE 4

Changing the duration of the synaptic input can switch the stability of synchrony. (A) Short duration synaptic input (τ_syn = 1.0 ms). (A1) Synaptic activation waveform over one cycle. (A2) Directly measured PRC (g_syn = 0.05 mS/cm²). The inset (more ...)

FIGURE 5

Stability of synchrony can also be switched in Wang and Buzsaki (WB) model neuron with thin action potentials. (A) Short duration synaptic input (τ_syn = 0.01 ms). (A1) Synaptic activation waveform over one cycle. (A2) Directly measured PRC (g (more ...)

We confirm the generality of these results with another example in Fig. 5, again using Wang and Buszaki neurons receiving inhibition for a brief and a longer duration synaptic time course. The rate of rise of synaptic activation was fast and constant. For the brief perturbation, the time constant τ_syn for the decay of the activation of the synaptic conductance was made arbitrarily fast (0.01 ms) in order to allow the width of the presynaptic action potential (identical here to that of the postsynaptic one) to control the duration of the elevation in activation of the synaptic conductance. The conductance is activated only while the presynaptic neuron is above threshold, effectively producing a square pulse with the width of an action potential for the sake of illustration. The brief perturbation (Fig. 5A1) is able to sample the initial anomalous region of positive slope in z_g(ϕ) (Fig. 5B) because it is of the same order as the action potential width (Fig. 5C). The thin action potentials cause the initial positive slope in the PRC (Fig. 5A2) to be barely discernable compared to the ML example in Fig. 4. However, this initial positive slope prevents two reciprocally coupled identical WB model neurons from synchronizing. For a slower decay time constant (τ_syn = 1.0 ms), the synaptic perturbation is broader than the action potential (Fig. 5E1) and is not able to preferentially sample the region of anomalous initial positive slope in z_g(ϕ). As a result, no initial anomalous region of positive slope is observed in the PRC (Fig. 5E2). The initial negative slope allows two reciprocally coupled identical WB model neurons in to synchronize in this case (Fig. 5E3). Hence the results are qualitatively similar for both the ML and WB models.

DISCUSSION

Two Equivalent Approaches to Weak Coupling

The theoretical framework for understanding phase resetting and its relationship to synchronization tendencies can be formulated in two ways, which we have shown in this manuscript to be formally equivalent. The approach based on the infinitesimal response with respect to current (z_i(ϕ)) is preferred by mathematicians because z_i(ϕ) is a scaled version of the solution to the adjoint linear system along the limit cycle. The adjoint can easily be calculated for a model system. Moreover, it does not have to be recalculated for different types of perturbations, such as those due to electrical synapses and the different types of excitatory and inhibitory chemical synapses. On the other hand, the approach based on infinitesimal response with respect to conductance (z_g(ϕ)) separates the perturbation from the system being perturbed as is usually assumed in systems theory (see Appendix (D)). For chemical synapses, different reversal potentials produce a different z_g(ϕ) which is useful for intuiting how the time course of the synapse will sample z(ϕ) to produce the PRC.

In order to use z_i(ϕ), we must either 1) obtain the perturbation in terms of the synaptic current i(ϕ)=g(ϕ)(V_post(ϕ)−V_syn) or 2) redefine the impulse response at a phase of ϕ by first scaling it by (V_post(ϕ)−V_syn). If we use the first approach, i(ϕ) as the input x(ϕ), we note that the input is not independent of the system, but requires a priori knowledge of one of the state variables of the system contained in the “black box” (see Fig.6 in Appendix (D)), V_post(ϕ). In weak coupling theory, the unperturbed value of this variable during a free running oscillation is used to assume a value for V(ϕ) at each phase. Similarly, if we use the second approach, a priori knowledge of one of the state variables of the system contained in the “black box” (see Fig.6 in Appendix (D)), V_post(ϕ) is required to calculate the impulse response at each phase. One motivation for this paper was to utilize a form for both the input and the impulse responses that, consistent with the systems theory approach, do not require direct measurement of any state variables within the “black box” of the system to compute the phase resetting curve, defined as the output given the input x(ϕ) and the infinitesimal phase varying response z(ϕ). This is accomplished by computing the response to an impulse in conductance directly, which produces z_g(ϕ), the infinitesimal phase resetting curve with respect to conductance instead of z_i(ϕ). The added advantage of using z_g(ϕ) is that it provides a better intuition with respect to what the effects of manipulations such as changing the frequency or the synaptic time course will be.

Correlation between PRC type and Excitability based on Bifurcation Structure

The discrimination between Type I and Type II excitability based upon the PRC is not as unambiguous as is implied by the statement that Type I PRCs consist only of either advances or delays whereas Type II consist of both (Ermentrout 1996). The emergence of spiking via a saddle node bifurcation on an invariant circle (SNIC) is associated with Type I excitability and with PRCs of predominantly one sign or the other (Ermentrout 1996); however, these PRCs will generally have an anomalous region (Oprisan and Canavier, 2002) if the perturbation is sufficiently brief to sample the anomalous region that always exists in the iPRC. Rinzel and Ermentrout (1998) stated that, with respect to applying an excitation to a neuron with Type I excitability, “In only a very small interval of time can the phase be delayed, and this is a general property of membranes that become oscillatory through a saddle node bifurcation”. This small interval corresponds to the action potential. We have shown that if parameters are adjusted so the action potential of a model with Type I excitability occupies a significant region of the cycle, its PRC becomes indistinguishable from a Type II PRC (Fig. 2B3). The assumption of thin action potentials should be made explicit when equating Type I excitability with Type I phase resetting. In order to avoid confusion, we should note that zero phase is frequently chosen in theoretical studies at the peak of the action potential (Ermentrout 1996) rather than at the threshold for synaptic activation as we have done here and is typically performed in experimental studies (Mancilla et al., 2007; Perkel et al., 1964; Preyer and Butera, 2005; Tateno and Robinson, 2007). This can shift the exact location of the anomalous region in the PRC.

With respect to spiking neurons, it is likely that the action potential occupies such a small fraction of the cycle that this anomalous region can often be ignored. The width of an action potential in the central nervous system (CNS) can vary from 180 μs to 4 ms (Bean, 2007). The time constants for fast EPSCs can range from the sub-millisecond range in specialized auditory nuclei to 2–5 msec in cortex (Destexhe, Mainen and Sejnowski, 1998), and fast IPSCs decay with a time constant of the order of 2 ms (Bartos et al., 2001). Thus it is conceivable that some PSCs in the CNS may have duration on the order of an action potential, making the findings illustrated in Figs. 2 and and33 more relevant to CNS neurons. For example, the WB model neuron exhibits thin action potentials in a Type I excitability regime using the parameters given in Appendix (B), but nonetheless can be made to exhibit a weakly Type II PRC that changes the observed synchronization properties. A more pronounced Type II PRC for oscillators exhibiting Type I excitability may also be observed in other contexts, such as for cardiac pacemakers, bursting neurons, or oscillators that subserve circadian rhythms. The Morris Lecar model with parameters as in Fig. 2 has relatively broad action potentials, and might be more analogous to the slow envelope oscillation that underlies bursting oscillations, for example. The PRCs of bursting neurons are more dependent upon the envelope rather than upon individual spikes (Oprisan et al., 2004; Sieling et al., 2009). The bifurcation structure that is most relevant to these PRCs should then be the one via which the slow envelope arises, not the bifurcations through which spiking is initiated and terminate (Izhikevich, 2000). If the oscillation that mediates the bursting envelope arises via a SNIC bifurcation, the resultant PRC could easily have prominent regions of both advances and delays if the hyperpolarizing phase of the envelope occupies a significant phase of the cycle. The fraction of the cycle occupied by the hyperpolarizing phase is frequency dependent, thus the degree to which these anomalous effects are manifested should be frequency dependent.

Experimentally, the distinction between Type I and Type II PRCs is also more complex than implied by the statement that Type I PRCs consist only of either advances or delays. For example, Tateno and Robinson (2007) characterized PRCs based on the absolute value of ratio of the peak advance to the peak delay or the ratio of the peak delay to the peak advance, whichever was smaller. If this ratio exceeded 0.175, the PRC was termed biphasic (Type II) otherwise it was termed monophasic (Type I). This type of fuzzy classification makes sense in view of the possibility that a very small anomalous region does not imply a Type II bifurcation structure. Furthermore, Tateno and Robinson (2007) found that whether the PRC was classified as monophasic or biphasic did not always strictly correspond with excitability type as measured by the frequency current curve. The presence of both advances and delays in the PRCs shown in Figs.2B3 and 3B3 does not imply that spiking arises via a Hopf bifurcation, as spiking clearly arises via a saddle node bifurcation in both cases. Instead, the change in shape of the PRC as the applied current becomes more depolarizing appears to reflect a transition from a Type I PRC shape to a Type II PRC shape as the Hopf bifurcation that terminates spiking is approached.

Synchronization Tendencies: Role of Synaptic Time Course

Generalizations with respect to the implications of the time course of the synaptic coupling on synchronization are frequently made. van Vreeswijk et al. (1994) stated that the infinitesimal PRC contained an anomalous region at only at phases corresponding to the action potential in which excitation produced delays and inhibition produced advances, which is a general characteristic of Type I excitability as shown in Figs. 4 and and5.5. They also stated in words the central idea that the convolution integral would produce under these circumstances an unstable slope in the PRC for synaptic interactions confined to the action potential and a stable one if the synaptic interactions sampled the region outside the action potential more heavily than the action potential. Our contribution is that we have made these calculations explicitly for a particular case in which the general predictions of van Vreeswijk et al. (1994) hold, that is, in which slow but not fast inhibition is synchronizing. Conversely, with respect to excitation, van Vreeswijk at al. (1994) stated that “In all the models we consider, perfectly synchronized firing cannot be produced by excitatory synaptic coupling unless the synaptic rise time is extremely rapid”. Their use of the term synaptic rise was linked to their use of an alpha coupling function, therefore instead of a rapid synaptic rise time, it is more accurate to say that, for a neuron with the appropriately shaped infinitesimal PRC, for excitation, synaptic durations shorter than an action potential favor synchrony. With respect to infinitesimal PRCs with different shapes, Hansel et al. (1995) point out that for Type II excitability, excitations applied just after the refractory period cause delays rather than advances, which can be synchronizing for Type II neurons when the interactions are “fast”, though presumably not so fast as to be confined to the duration of an action potential. The theoretical framework in this paper makes the role of the speed of interactions more precise and generalizable.

Another implication of this work is that the same input may produce different synchronization tendencies in different frequency ranges; an input that is brief at low frequency is not so brief relative to the cycle period at higher frequencies.

Summary

In this paper, we presented a theoretical framework that clarifies some misconceptions in the literature such as the idea that Type I excitability/bifurcation type always implies a monophasic (Type I) phase resetting curve (Ermentrout 1996) and that synaptic rise time determines synchronization properties (van Vreeswijk et al., 1994). We showed how changing the frequency of a neural oscillator can switch the observed PRC type from monophasic (Type I) to biphasic (Type II) simply by increasing the fraction of the cycle occupied by the action potential as shown in Figs. 2 and and3.3. This is important because monophasic resetting for Type I excitability is only guaranteed at phases that fall in the interspike interval, and not during the action potential. We also showed how changing the duration of a perturbation can change the slope of PRC at the locking point and thereby switch the stability of the synchronous state as shown in Figs. 4 and and5.5. The real determinant of whether excitation or inhibition stabilizes synchrony is the shape of the PRC, and this shape is strongly influenced by both the intrinsic properties of the neuron and the properties of the synapse.

Ackowledgements

This work was supported by NIH NS54281 under the CRCNS program to CCC. SA would like to thank Dr. Selva K. Maran and Fred Sieling for valuable discussions.

Appendix (A): Morris Lecar Model Neuron Equations and Parameters

The current balance equation for the Morris-Lecar model (ML model) neuron (Morris and Lecar, 1981; Rinzel and Ermentrout, 1998) is:

where the membrane capacitance C_m = 1μF/cm², V_m is the cell membrane voltage in millivolts and t is time in milliseconds. The leak current is given by I_L =g_L(V_m−E_L). The calcium current is given by I_Ca =g_Cam_∞(V_m)(V_m −V_Ca). The steady state activation is :

where V₁= −1.0 mV and V₂=15.0 mV. The potassium current is given by I_K =g_Kw(V_m−V_K). The rate equation for the activation variable w in the expression for the potassium current is:

where

was set to 0.2 with the steady state activation amplitude as:

and activation rate as:

with V₃=10 mV and V₄=14.5 mV. The reversal potentials V_Ca, V_K and V_L were set to 100, −70.0 and −50 mV, respectively. The maximal potassium (g_K) and leak (g_L) conductances were set to 2.0 and 0.5 mS/cm², respectively. For calcium, the maximal conductance maximal (g_Ca) was set to 1.0 mS/cm² and I_stim was set at 9.0 μA/cm². Unless otherwise stated, the values for the various parameters for the Type I excitability regime were equal to those given above and were taken from Ermentrout, 2002. The synaptic current is given by I_syn =g_syns(V_m −V_syn), where g_syn is the maximum synaptic conductance and V_syn is equal to −75.0 mV for inhibitory synaptic connectivity and equal to 0.0 mV for excitatory synaptic connectivity. The rate of change of the gating variable s is given by ds/dt = α T(V_pre)(1−s) − s/τ_syn with T(V_pre)=1/[1+exp(−(V_pre −V_half)/2)], where V_pre is the voltage of the presynaptic cell, V_half is the synaptic half activation voltage (was taken to be 0mV if the threshold voltage was set to −14.0 mV and was taken to be 28mV if the threshold voltage was set to 14.0 mV), α=6.25ms⁻¹ is the rate constant of the synaptic activation (Bartos et al. 2001), and τ_syn is the synaptic decay time constant and was set to 1.0 ms or 3.0 ms depending on brief or long lasting synapese, respectively.

Appendix (B): Wang and Buzsaki Model Neuron Equations and Parameters

The current balance equation for the Wang and Buzsaki (WB) model neuron (Wang and Buzsaki, 1996) is:

where the capacitance C = 1μF/cm², V is the cell membrane voltage in millivolts and t is time in milliseconds. The leak current is given by I_L =g_L(V−V_L). The sodium current is given by I_Na=g_Nam_∞³h(V−V_Na). The steady state activation m_∞=α_m/(α_m+β_m) where α_m(V)= −0.1(V+35)/{exp[−0.1(V+35)]−1}and β_m(V)=4exp[−(V+60)/18]. The rate equation for the inactivation variable h in the expression for sodium current is:

where

=5. (The symbol was changed from ϕ to [var phi]

to avoid confusion with the symbol for phase). The rate constants for the inactivation variable h are given by α_h(V)= 0.07exp[−(V+58)/20] and β_h(V)= 1/{exp[−0.1(V+28)]+1}. The potassium current is given by I_k =g_kn⁴(V−V_K), where the activation variable n satisfies the following equation:

where the rate constants for n are α_n(V)= −0.01(V+34)/{exp[−0.1(V+34)]−1} and β_n(V)=0.125exp[−(V+44)/80]. The reversal potentials V_Na, V_K and V_L were set to 55, −90 and −65 mV, respectively. The maximal sodium (g_Na), potassium (g_K) and leak (g_L) conductances were set to 35, 9 and 0.1 mS/cm², respectively. I_stim is the applied current and was set at 0.5 μA/cm². Unless otherwise stated, the values for the various parameters were equal to those given above. The synaptic coupling current was the same as in the case of the ML model neuron with threshold voltage set to −14.0 mV and V_half equal to 0.0 mV. The synaptic decay time constant τ_syn was set to 1.0 ms for the case of long duration input (see Fig.6E1) and 0.01 ms for short duration input (see Fig. 6A1).

Appendix (C): Reinterpretation of a Phase Oscillator as a Linear Time Varying System

A system responds to a signal (the input) by producing another signal (the output). For a linear system (Oppenheim et al., 1983), the output, y(t), to any arbitrary signal, x(t), can be predicted using only the measured response h_τ(t) of the system to a single unit impulse at time `τ'. This impulse is designated as the unit delta function δ(t−τ) for continuous systems, which is zero everywhere except at time τ but has unit area. With no simplifications, y(t)=∫x(τ)h_τ(t)dτ. For time invariant systems, h_τ(t) = h₀(t−τ) = h(t−τ), i.e. the response to an unit impulse applied at time τ is identical to the response at any other time except but shifted in time (see Fig. 6A). Therefore for a linear time invariant system the general expression simplifies to a convolution integral y(t)=∫x(τ) h(t−τ)dτ. If x(t) is represented as train of scaled and shifted unit impulses such that x(t)=∫x(τ) δ(t−τ) dτ, then the superposition of the impulse responses to each impulse in the train gives the output by the assumption of linearity.

The average phase difference theory for weakly coupled oscillators (Ermentrout and Kopell, 1991) stands alone as originally formulated, but can also be understood in terms of linear systems theory. First, we assume a periodic oscillator that repeats the same closed path in the state space of the system. This closed path is called a stable limit cycle, and all nearby trajectories are attracted to the limit cycle. The independent variable is phase ϕ, which is defined in this context as position along the limit cycle modulo 1. The output of the system is the phase resetting curve H(ϕ), which is not a time series, but instead is defined at each point along the limit cycle. Since the state of an oscillator is presumed to be completely characterized by the phase, such an oscillator is called a phase oscillator. The phase σ at which an impulse perturbation δ(ϕ−σ) is applied changes the impulse response h_σ(ϕ−σ), therefore phase oscillators are not invariant with respect to the independent variable. The second assumption is that the system is only weakly nonlinear such that the principle of superposition still applies: the output produced by any two inputs is equal to the sum of the outputs produced in response to each input applied separately. A third assumption is that the phase resetting is instantaneous, such that the response to an impulse δ(ϕ−σ) at a given phase σ is itself an impulse. The amplitude of this impulse response depends on the phase at which the input impulse is applied: h_σ(ϕ) = ∫z(ϕ)δ(ϕ−σ)dϕ = z(σ) (see Fig. 6B). z(ϕ) is the infinitesimal phase response curve. An arbitrary input perturbation x(ϕ) can also be defined as a function of position on the limit cycle over the interval 0 to 1 so that the independent variable is consistent across the system. In order to obtain the phase response curve at each position σ along the limit cycle, the input x(ϕ) is shifted with respect to the position on the limit cycle so that it becomes x(ϕ+σ) for an input that is applied beginning at a phase of σ. The fourth and final assumption is that the phase need not be updated as a result of the weak phase resetting that occurs earlier in the cycle, such that the total phase resetting due to a weak input x(ϕ) applied at a phase of σ is H(σ) = ∫ x(ϕ+σ)z(ϕ)dϕ evaluated using the unperturbed phase over a cycle period from ϕ=0 to ϕ=1, which is a circular convolution of periodic signals. This can equivalently be written as H(σ) = ∫ x(ϕ)z(ϕ−σ)dϕ for comparison with the time-invariant linear system. To summarize, for linear time invariant system, a convolution results because the impulse response is invariant, whereas for a phase oscillator, a convolution results because the impulse response is instantaneous.

Appendix (D): Sources of Error in the Weak Coupling Estimation of the Directly Measured PRC

To compute the PRC estimated by the iPRC approach with respect to conductance pulse (z_g(ϕ)) which is sampled at a finite number (n) of phases (see Methods), the following convolution sum for discrete systems is used.

In the above eqn., Δt represents the interval between the pulses of the input train, g_syn(ϕ−σ[n]). Δσ[n] refers to the phase resetting due to each pulse of the train at phase σ[n]. In order to illustrate the error introduced by not updating the phase due to phase resetting earlier in the cycle, we also used the following form for comparison in which the phase resetting due to previous pulses in the train is incorporated in the calculation of the phase resetting due to subsequent pulses:

Figure 7A illustrates that for a conductance strength of 0.005 mS /cm², which is five times that used in the example given in Fig. 1B, the first order PRC computed by the direct method (see Methods) is not an exact match to H_g(ϕ). On the other hand, Ĥ_g(ϕ), which updates the phase due to resetting earlier in the cycle, is a better match to the directly measured PRC. Achuthan and Canavier (2009) showed that as the strength of the resetting is increased, the principle of superposition breaks down in that resetting no longer scales with the conductance. Fig. 7A shows that the assumption that the phase need not be updated to account for resetting earlier in the cycle breaks down before the superposition assumption, at least in this example. The remaining error is quite small and occurs only at late phases. Since the circular sum in both H_g(ϕ) and Ĥ_g(ϕ) is performed over one full cycle period starting at the phase at which the perturbation is applied, the length of both the cycle that contains the start of the perturbation and that of the subsequent cycle can be affected. Therefore, the result of the convolution actually approximates the sum of the first and second order resetting, rather than the first order resetting alone. According, the remaining error can be removed by comparing Ĥ_g(ϕ) (see Fig.7B) to the sum of the first and second order resetting obtained by the direct method, instead of to the first order resetting alone.

FIGURE 7

Sources of error in the weak coupling estimation of the directly measured PRC. (A) Error introduced by not updating the phase continuously. The first order PRC₁(solid black curve) obtained by the direct method matches Ĥ_g(ϕ) (black dashed (more ...)

Footnotes

¹Fig 6A of Maran and Canavier mistakenly implied that the WB model had a significant instantaneous second order reset. The pulse height but not the width went to zero in that figure; when the width is also made to go to zero, the second order resetting essentially disappears.

²Izhikevich (2007) pointed out that Hodgkin's definition of Type I excitability is only strictly satisfied by a saddle node bifurcation on an invariant limit cycle in which the saddle node collides with the quiescent stable fixed point and the limit cycle simultaneously.

REFERENCES

Achuthan S, Canavier CC. Phase resetting curves determine synchronization, phase locking, and clustering in networks of neural oscillators. J. Neurosci. 2009;29(16):5218–5233. [PMC free article] [PubMed]
Acker CD, Kopell N, White JA. Synchronization of strongly coupled excitatory neurons: relating network behavior to biophysics. J. Comput. Neurosci. 2003;15(1):71–90. [PubMed]
Bartos M, Vida I, Frotscher M, Geiger JRP, Jonas P. Rapid signaling at inhibitory synapses in a dentate gyrus interneuron network. J. Neurosci. 2001;21(8):2687–2698. [PubMed]
Bean BP. The action potential in mammalian central neurons. Nature Rev. Neurosci. 2007;8:451–465. [PubMed]
Brown E, Moehlis J, Holmes P. On the phase reduction and response dynamics of neural oscillator populations. Neural Comput. 2004;16(4):673–715. [PubMed]
Buzsaki G. Rhythms of the brain. Oxford University Press Inc.; New York: 2006.
Canavier CC, Butera RJ, Dror RO, Baxter DA, Clark JW, Byrne JH. Phase response characteristics of model neurons determine which patterns are expressed in a ring circuit model of gait generation. Biol. Cyber. 1997;77(6):367–380. [PubMed]
Canavier CC, Baxter DA, Clark JW, Byrne JH. Control of multistability in ring circuits of oscillators. Biol. Cyber. 1999;80:87–102. [PubMed]
Destexhe A, Mainen ZF, Sejnowski TJ. Kinetic models of synaptic transmission. In: Koch C, Segev I, editors. Methods in neuronal modeling from ions to networks. The MIT Press; Cambridge, MA: 1998.
Dror RO, Canavier CC, Butera RJ, Clark JW, Byrne JH. A mathematical criterion based on phase response curves for stability in a ring of coupled oscillators. Biol.Cybern. 1999;80:11–23. [PubMed]
Ermentrout GB, Kopell N. Frequency plateaus in a chain of weakly coupled oscillators, I. SIAM J. Math. Anal. 1984;15(2):215–237.
Ermentrout GB, Kopell N. Oscillator death in systems of coupled neural oscillators. SIAM J. Appl. Math. 1990;50(1):125–146.
Ermentrout GB, Kopell N. Multiple pulse interactions and averaging in coupled neural oscillators. J. Math Biol. 1991;29(3):195–217.
Ermentrout GB. Type I membranes, phase resetting curves and synchrony. Neural Comput. 1996;8(5):979–1001. [PubMed]
Ermentrout GB. Simulating, analyzing, and animating dynamical systems: A guide to XPPAUT for researchers and students. SIAM; Philadelphia: 2002.
Fell J, Klaver P, Lehnertz K, Grunwald T, Schaller C, Elger CE, Fernandez G. Human memory formation is accompanied by rhinal-hippocampal coupling and decoupling. Nat. Neurosci. 2001;4(12):1259–1264. [PubMed]
Galan RF, Ermentrout GB, Urban NN. Efficient estimation of phase resetting curves in real neurons and its significance for neural modeling. Phys. Rev. Lett. 2005;94(15):158101. [PMC free article] [PubMed]
Goel P, Ermentrout GB. Synchrony, stability, and firing patterns in pulse-coupled oscillators. Physica D. 2002;163(3– 4):191–216.
Guckenheimer J, Holmes P. Nonlinear oscillations, dynamical systems and bifurcations of vector fields. Springer-Verlag, Inc.; New York: 1983.
Gutkin BS, Ermentrout GB, Reyes AD. Phase-response curves give the responses of neurons to transient inputs. J. Neurophysiol. 2005;94:1623–1635. [PubMed]
Hammond C, Bergman H, Brown P. Pathological synchronization in Parkinson's disease: networks, models and treatments. Trends Neurosci. 2007;30(7):357–364. [PubMed]
Hansel D, Mato G, Meunier C. Synchrony in excitatory neural networks. Neural Comput. 1995;7(2):307–337. [PubMed]
Hodgkin AL. The local electric changes associated with repetitive action in a non-medullated axon. J. Physiol. 1948;107(2):165–181. [PubMed]
Huguenard JR, McCormick DA. Thalamic synchrony and dynamic regulation of global forebrain oscillations. Trends Neurosci. 2007;30(7):350–356. [PubMed]
Hutchison WD, Dostrovsky JO, Walters JR, Courtemanche R, Boraud T, Goldberg J, Brown P. Neuronal oscillations in the basal ganglia and movement disorders: Evidence from whole animal and human recordings. J. Neurosci. 2004;24(42):9240–9243. [PubMed]
Hoppensteadt FC, Izhikevich EM. Weakly connected neural networks. Springer-Verlag; New York: 1997.
Izhikevich EM. Neural excitability, spiking and bursting. Int. J. of Bifurcat. and Chaos. 2000;10:1171–1266.
Izhikevich EM. Dynamical systems in neuroscience: the geometry of excitability and bursting. The MIT Press; Cambridge, MA: 2007.
Khoo MCK. Physiological control systems: analysis, simulation and estimation. IEEE Press; New York, NY: 2000.
Mancilla JG, Lewis TJ, Pinto DJ, Rinzel J, Connors BW. Synchronization of electrically coupled pairs of inhibitory interneurons in neocortex. J. Neurosci. 2007;27(8):2058–2073. [PubMed]
Maran SK, Canavier CC. Using phase resetting to predict 1:1 and 2:2 locking in two neuron networks in which firing order is not always preserved. J. Comput. Neurosci. 2008;24(1):37–55. [PMC free article] [PubMed]
Morris C, Lecar H. Voltage oscillations in the barnacle giant muscle fiber. Biophys. J. 1981;35(1):193–213. [PubMed]
Netoff TI, Acker CD, Bettencourt JC, White JA. Beyond two-cell networks: experimental measurement of neuronal responses to multiple synaptic inputs. J. Comput. Neurosci. 2005a;18(3):287–295. [PubMed]
Netoff TI, Banks MI, Dorval AD, Acker, C. D, Haas JS, Kopell N, White JA. Synchronization in hybrid neuronal networks of the hippocampal formation. J. Neurophysiol. 2005b;93:1197–1208. [PubMed]
Oppenheim AV, Willsky AS, Young IT. Signals and systems. Prentice-Hall, Inc.; Englewood Cliffs, NJ: 1983.
Oprisan SA, Canavier CC. The influence of limit cycle topology on the phase resetting curve. Neural Comput. 2002;14(5):1027–1057. [PubMed]
Oprisan SA, Prinz AA, Canavier CC. Phase resetting and phase locking in hybrid circuits of one model and one biological neuron. Biophys. J. 2004;87(4):2283–2298. [PubMed]
Perkel DH, Schulman JH, Bullock TH, Moore GP, Segundo JP. Pacemaker neurons: effects of regularly spaced synaptic input. Science. 1964;145(3627):61–63. [PubMed]
Preyer AJ, Butera RJ. Neuronal oscillators in aplysia californica that demonstrate weak coupling in vitro. Phys. Rev. Lett. 2005;95(13):138103. [PubMed]
Rinzel J, Ermentrout GB. Analysis of neural excitability and oscillations. In: Koch C, Segev I, editors. Methods in neuronal modeling from ions to networks. The MIT Press; Cambridge, MA: 1998.
Rodriguez E, George N, Lachaux JP, Martinerie J, Renault B, Varela FJ. Perception's shadow: long-distance synchronization of human brain activity. Nature. 1999;397(6718):430–433. [PubMed]
Sharp AA, O'Neil MB, Abbott LF, Marder E. The dynamic clamp - artificial conductances in biological neurons. Trends Neurosci. 1993a;16(10):389–394. [PubMed]
Sharp AA, O'Neil MB, Abbott LF, Marder E. Dynamic clamp: computer-generated conductances in real neurons. J. Neurophysiol. 1993b;69:992–995. [PubMed]
Sieling FH, Canavier CC, Prinz AA. Predictions of phase-locking in excitatory hybrid networks: excitation does not promote phase-locking in pattern-generating networks as reliably as inhibition. J. Neurophysiol. 2009;102:69–84. [PubMed]
Stelt O, van der, Belger A, Lieberman JA. Macroscopic fast neuronal oscillations and synchrony in schizophrenia. Proc. Natl. Acad. Sci. USA. 2004;101(51):17567–17568. [PubMed]
Tateno T, Robinson HPC. Phase resetting curves and oscillatory stability in interneurons of rat somatosensory cortex. Biophys. J. 2007;92(2):683–695. [PubMed]
Timofeev I, Steriade M. Neocortical seizures: initiation, development and cessation. Neuroscience. 2004;123(2):299–336. [PubMed]
Traub RD, Jefferys JG. Are there unifying principles underlying the generation of epileptic afterdischarges in vitro? Prog. Brain Res. 1994;102:383–394. [PubMed]
Uhlhaas PJ, Singer W. Neural synchrony in brain disorders: relevance for cognitive dysfunctions and pathophysiology. Neuron. 2006;52(1):155–168. [PubMed]
Van Vreeswijk C, Abbott LF, Ermentrout GB. When inhibition not excitation synchronizes neural firing. J. Comput. Neurosci. 1994;1(4):313–321. [PubMed]
Wang XJ, Buzsaki G. Gamma oscillation by synaptic inhibition in a hippocampal interneuronal network model. J. Neurosci. 1996;16(20):6402–6413. [PubMed]
Winfree AT. The geometry of biological time. Springer; New York: 2001.