## Abstract

One effect of any external perturbations, such as presynaptic inputs, received by limit cycle oscillators when they are part of larger neural networks is a transient change in their firing rate, or phase resetting. A brief external perturbation moves the figurative point outside the limit cycle, a geometric perturbation that we mapped into a transient change in the firing rate, or a temporal phase resetting. In order to gain a better qualitative understanding of the link between the geometry of the limit cycle and the phase resetting curve (PRC), we used a moving reference frame with one axis tangent and the others normal to the limit cycle. We found that the stability coefficients associated with the unperturbed limit cycle provided good quantitative predictions of both the tangent and the normal geometric displacements induced by external perturbations. A geometric-to-temporal mapping allowed us to correctly predict the PRC while preserving the intuitive nature of this geometric approach.

## 1.  Introduction

Since the first recordings from Aplysia neuron R15 (Alving, 1968), we have learned that some neurons are intrinsic oscillators. Neural oscillators that produce periodic rhythms can be modeled as limit cycle oscillators (Baxter et al., 2000). Networks of nonlinear oscillators have attracted a great deal of research due to their theoretical and practical importance (Ermentrout & Kopell, 1985, 1986; Ermentrout, 1994; Izhikevich, 2000). When acting as a component of a neural network, nonlinear oscillators can be characterized by their phase resetting curve (PRC) (Abramovich-Sivan & Akselrod, 1998a, 1998b; Dror, Canavier, Butera, & Byrne, 1999; Ermentrout, 1996; Pinsker, 1977). The PRC measures the transient change in the period of the limit cycle when the same perturbation is applied at different phases during the cycle. Depending on the strength and the type of coupling (inhibitory or excitatory) among neural oscillators, as well as the topology of the neural network (e.g., ring, all-to-all), neural oscillators can develop stationary phasic relationships throughout the entire network called phase-locked modes. The PRC is used for predicting phase-locked modes in neural networks even when the explicit mathematical models of component neurons are not available (Oprisan, Prinz, & Canavier, 2004). In such cases, the PRC simply tabulates the experimentally measured transient changes in the firing rate of an isolated neuron (open loop) subject to an external stimulus identical to that received by the same neuron as part of a fully connected network (closed loop) (Oprisan & Boutan, 2008; Oprisan, 2009). Memory storage (Haenschel, Linden, Bittner, Singer, & Hanslmayr, 2010; Rutishauser, Ross, Mamelak, & Schuman, 2010) and information processing (Zacksenhouse & Ahissar, 2006; Cheron et al., 2007; Wingerden, Vinck, Lankelma, & Pennartz, 2010) are important examples of phase-locking processes that take place in neural networks.

The term phase has been used in the literature with different meanings. Winfree (1980) defined the phase, , as the elapsed time, t, measured from an arbitrary reference divided by the intrinsic period of oscillations, Pi: . We called the temporal phase to distinguish it from the geometrical phase, (see Oprisan & Canavier, 2002, 2006; Oprisan, 2012a, 2012b). The geometrical phase, , is the normalized distance measured along the limit cycle, that is, phase space distance divided by the total phase space length of the limit cycle. The Euclidean distance between two points along is given by
1.1
where is the length of covered by the figurative point over a period of oscillation. The geometrical phase is a monotonically increasing function of the elapsed time t, or the equivalent temporal phase (see Figure 1B). The relationship is an invertible map between the temporal and the geometrical phases, which can be used to estimate the PRC based on the numerically computed geometric perturbations of the limit cycle (see Oprisan & Canavier, 2002, for details). Alternatively, Murray (1993) and Hoppensteadt and Izhikevich (1997) defined the angular phase as the angle that parameterizes the unit circle that maps any limit cycle oscillator. For a constant angular velocity on a circular limit cycle, the temporal and angular phases are equivalent.
Figure 1:

Temporal and geometric phase of oscillators. (A) A typical limit cycle of a Bonhoeffer–van der Pol equation (BvP) neural oscillator (solid line) with phases marked every 5% of the intrinsic period (filled circles). The fast variable nullcline, f1=0, (dashed line) and the slow variable nullcline f2=0 can have multiple crossings called equilibrium points, steady states, or fixed points. (B) The slope of the geometric distance, , traveled by the figurative point along the limit cycle is proportional to the local flow. A current/conductance perturbation pushes the figurative point outside the (unperturbed) limit cycle, which subsequently returns to at a geometrical phase difference (vertical arrow) from its twin figurative point that traveled only along the unperturbed limit cycle . The geometric phase resetting is mapped into a temporal phase difference (horizontal arrow) by the invertible map . (C) Numerically computed PRCs obtained by directly integrating the model's equations (solid line) and the predicted PRCs based on the invertible geometric-to-temporal phase maps (dashed line) for a of the BvP model.

Figure 1:

Temporal and geometric phase of oscillators. (A) A typical limit cycle of a Bonhoeffer–van der Pol equation (BvP) neural oscillator (solid line) with phases marked every 5% of the intrinsic period (filled circles). The fast variable nullcline, f1=0, (dashed line) and the slow variable nullcline f2=0 can have multiple crossings called equilibrium points, steady states, or fixed points. (B) The slope of the geometric distance, , traveled by the figurative point along the limit cycle is proportional to the local flow. A current/conductance perturbation pushes the figurative point outside the (unperturbed) limit cycle, which subsequently returns to at a geometrical phase difference (vertical arrow) from its twin figurative point that traveled only along the unperturbed limit cycle . The geometric phase resetting is mapped into a temporal phase difference (horizontal arrow) by the invertible map . (C) Numerically computed PRCs obtained by directly integrating the model's equations (solid line) and the predicted PRCs based on the invertible geometric-to-temporal phase maps (dashed line) for a of the BvP model.

Neural oscillators are dynamical systems described by a nonlinear flow F:
1.2
where defines a flow in and is a control parameter in —for example, bias current, ionic channel conductances, or half activation voltage of ionic channels. In all neurobiological models, x1 designates the membrane potential, and the other variables are usually ascribed to gating variables and ionic concentrations, for example (Brown, Moehlis, & Holmes, 2004; Hoppensteadt & Izhikevich, 1997; Izhikevich, 2007). A stable limit cycle, , of the dynamical system given by equation 1.2 is a periodic solution x(t)=x(t+Pi) that attracts all points in a neighborhood of as .

As a concrete example, we show the phase-plane (x1, x2) orbits of a two-dimensional Bonhoeffer–van der Pol equations (BvP) (FitzHugh, 1961) (see equation 3.1 and section 3 for details regarding the model equations) with a solid line in Figure 1A. Equally spaced temporal marks every 5% of Pi are shown in Figure 1A with filled circles. The phase reference is arbitrary, and in Figure 1A, it corresponds to crossing of the minimum membrane voltage with positive slope.

The geometric-to-temporal phase map, such as the one shown in Figure 1B for the BvP model, allowed us to estimate the PRC by converting the numerically computed geometric displacement that occurs between the figurative point that traveled on a perturbed trajectory outside the limit and the corresponding unperturbed figurative point to a temporal phase difference . We previously derived quantitative expressions that estimated the geometric phase difference due to a perturbation of the membrane voltage. It was previously established that the tangent geometrical displacement induced by a time-dependent injected current i(t) is given by , where Cm is the membrane capacitance, is the duration of the stimulus, f1 is the rate of change of the membrane potential, and ||F(t)|| is the Euclidean norm of the vector field (see equations 2.4 to 2.6 in Oprisan & Canavier, 2002, but also Oprisan & Canavier, 2006; Oprisan, 2012a, 2012b). For planar dynamical systems, the contribution of the normal geometrical displacement to the PRC was estimated by (see equation 3.3 in Oprisan & Canavier, 2002, but also Oprisan & Canavier, 2006; Oprisan, 2012a, 2012b). Using the numerically computed tangent and normal projections of the geometric displacement with the above formulas and the invertible map obtained from the BvP model's equations (see Figure 1B), we numerically estimated the temporal phase resetting (see Figure 1C). The geometric-to-temporal mapping method worked very well for type 1 excitable cells for which the periodic oscillations emerge via a saddle node of invariant circle (SNIC) bifurcation (see Ermentrout, 1996; Hoppensteadt & Izhikevich, 1997; Izhikevich, 2007, for detailed classifications of excitable cells according to their bifurcation structures). For type 1 excitable cells, the PRC is unimodal and corresponds to either only advances or only delays of the subsequent responses (see, e.g., Oprisan & Canavier, 2002). However, the same geometric-to-temporal mapping method was less successful when applied to type 2 excitable cells (see Figure 8 in Oprisan & Canavier, 2002, and Figure 1C in this letter) for which the periodic oscillations emerge via an Andronov-Hopf bifurcation (HB) (see Ermentrout, 1996; Hoppensteadt & Izhikevich, 1997; Izhikevich, 2007). For type 2 excitable cells, the PRC is bimodal (see Figure 1C for the BvP model and additional examples for the Morris-Lecar model neuron in Oprisan & Canavier, 2002). The reason that the geometric-to-temporal map worked so well for type 1 excitability class is that the phase-space position and the shape of the limit cycle of such oscillators are relatively insensitive to the bias current perturbations. Therefore, for a type 1 model neuron, a perturbation produces a negligible displacement of a figurative point along the direction normal to the limit cycle (i.e., ), since the two limit cycles (unperturbed and perturbed) overlap almost perfectly (see, e.g., Figure 5 in Oprisan & Canavier, 2002). As a result, in the case of type 1 neural oscillators, the geometrical phase difference is almost entirely due to sudden tangent displacements of the figurative point along the limit cycle. Previously derived quantitative formulas (see equations 2.4 to 2.6 in Oprisan & Canavier, 2002) captured very well the tangent geometric displacement and therefore produced accurate PRC predictions (see, e.g., Figure 5 in Oprisan & Canavier, 2002). On the other hand, for a type 2 neural oscillator, the phase-space position and the shape of the limit cycle change significantly with the bias current (see Figure 6 in Oprisan & Canavier, 2002, for a detailed example). Therefore, for a type 2 model neuron, a bias current perturbation produces a large displacement of the figurative point along the direction normal to the limit cycle since the two limit cycles (unperturbed and perturbed) are significantly shifted with respect to each other (see, e.g., Figure 6 in Oprisan & Canavier, 2002). As a result, in the case of type 2 neural oscillators, a better understanding of the normal component of the geometrical phase perturbation is crucial for a more accurate quantitative estimation of the PRC. Note that the geometric approach to the PRC prediction does not make any assumptions regarding the strength or the temporal evolution of the perturbation, and the method has been used for predicting the PRCs of neural oscillators subject to arbitrary time-dependent stimuli (Oprisan & Canavier, 2002, 2006; Maran, Sieling, Demla, Prinz, & Canavier, 2011). The geometrical-to-temporal phase mapping is invariant to changes of variables and can be obtained using either the model's equations (if known) or performing a time series delay embedding (see Oprisan, Thirumalai, & Canavier, 2003, for a detailed example pertaining to the stomatogastric ganglion recordings).

To conclude this section, we notice that the response of a dynamical system to external perturbations is conveniently measured in the temporal domain; for example, brief electric stimuli reset the temporal phase of a neural oscillator and produce a measurable change in the firing period of a neuron. The reality is that electrical stimuli do not directly interact with the abstract temporal dimension of a neural oscillator. The electric stimuli produce measurable changes in the phase-space trajectory of a neuron by directly changing (some of the) intrinsic variables, such as the fraction of ionic channels open. As a result, an external perturbation changes the local topology of the limit cycle of the neural oscillator and induces a total geometric phase resetting . The connection between the geometry of the phase space and the temporal dimension of a dynamical system is given by the geometric-to-temporal phase map (see Oprisan & Canavier, 2002, 2006), which is analyzed in this letter from a different perspective that relates and the stability coefficients of the limit cycle.

## 2.  Methods

Why is the PRC prediction for type 2 neurons based on the geometrical-to-temporal mapping approach (see the dashed line in Figure 1C) not as accurate as the prediction given by the same method for type 1 neurons? The obvious answer is that the estimation of the relaxation time due to the normal displacement given by the above formula is problematic (see also equation 3.3 in Oprisan & Canavier, 2002). We believe that a better qualitative understanding and a potentially significant quantitative improvement in the PRC prediction can be achieved by a more realistic estimation of the tangent and normal perturbations based on the stability coefficients associated with the limit cycle in the moving reference frame of the figurative point. Briefly, the method we use here considers that the figurative point on the limit cycle is perturbed by a vector (in the direction of membrane potential) and that the local eigenvalues on the limit cycle dictate different relaxation speeds in the directions normal and tangent to the limit cycle at the point of the perturbation. In order to compute the PRC, we determined the intersection of the relaxation trajectory with the unperturbed limit cycle and mapped the geometric distance between the return point and the unperturbed figurative point into a temporal phase resetting .

When investigating the motion of figurative points along a limit cycle , it is more convenient to use a moving reference frame instead of the fixed frame of equation 1.2. The goal of converting from the fixed coordinated x to the moving frame can be achieved by using a unitary transform , with , where the superscript T means transposed and vi are orthogonal vectors (Ali & Menzinger, 1999; Bellman, 1997; Shores, 2007). For a given N-dimensional flow F, the unitary matrix U is not uniquely defined. However, the tangent direction to the limit cycle can always be chosen as . The remaining N−1 normal vectors to the limit cycle cannot be uniquely defined, except for the fact that they must be orthogonal to each other.

Without reducing the generality of our results, we explained how the geometric approach to the PRC prediction works for planar vector flows. In this case, the tangent vector is , and there are two possible choices for the normal vector v2. Since we selected clockwise motion around the limit cycle (see Figure 2), the normal vector v2 is given by . As a result, the unitary matrix is
2.1
Figure 2:

Fixed versus moving reference frame in planar dynamical systems. The position of a figurative point with respect to the fixed reference frame (x1, x2) is given by the position vector r(x1, x2). The filled circles represent temporal markers 5% of the intrinsic period apart. The moving reference frame (dashed lines) is defined by the unit vectors (v1, v2) and has one axis tangent, xt, and the other normal, xn, to the limit cycle at any point. A brief perturbation (thick arrow) moves the figurative point A (filled square) outside the limit cycle to a point A’ (star). After the perturbation is over, it takes seconds for the perturbed trajectory starting at A’ (dashed-dotted line) to relax to an neighborhood of the unperturbed limit cycle (point B’). During the same time interval, , the figurative point A would have traveled along the unperturbed limit cycle (continuous line) to point B. The spatial, or geometric, distance between points B and B’ maps into a temporal phase difference and determines the PRC in response to the applied perturbation.

Figure 2:

Fixed versus moving reference frame in planar dynamical systems. The position of a figurative point with respect to the fixed reference frame (x1, x2) is given by the position vector r(x1, x2). The filled circles represent temporal markers 5% of the intrinsic period apart. The moving reference frame (dashed lines) is defined by the unit vectors (v1, v2) and has one axis tangent, xt, and the other normal, xn, to the limit cycle at any point. A brief perturbation (thick arrow) moves the figurative point A (filled square) outside the limit cycle to a point A’ (star). After the perturbation is over, it takes seconds for the perturbed trajectory starting at A’ (dashed-dotted line) to relax to an neighborhood of the unperturbed limit cycle (point B’). During the same time interval, , the figurative point A would have traveled along the unperturbed limit cycle (continuous line) to point B. The spatial, or geometric, distance between points B and B’ maps into a temporal phase difference and determines the PRC in response to the applied perturbation.

As discussed briefly in section 1 and extensively documented in Oprisan and Canavier (2002, 2006), for a type 1 neural oscillator, the most significant displacement of the figurative point due to a perturbation is along the direction tangent to the unperturbed limit cycle. Type 2 oscillators have a limit cycle that significantly shifts in the phase space as a result of external perturbations leading to large normal displacements of the figurative point relative to the unperturbed limit cycle. In general, perturbations of limit cycle oscillators lead to both normal and tangent displacements from the limit cycle (see the thick arrow in Figure 2). For stable limit cycles, such perturbations eventually die out and the figurative point returns to the limit cycle , or an -neighborhood of it, with an arbitrarily small positive number. The geometric distance, measured in the phase space of the dynamical system, between the figurative point that travels along the unperturbed limit cycle (from A to B in Figure 2) and the returning point after the perturbation dies out (point B′ in Figure 2) represents the amount of geometrical phase resetting induced by the perturbation (see Figure 2).

For small perturbations, , from the unperturbed limit cycle , equation 1.2 reduces to
2.2
where
2.3
is the Jacobian matrix evaluated along the limit cycle .
Using the unitary transform from the fixed reference frame x=(x1, x2)T to the moving reference frame , that is, and equation 1.2, we get (Ali & Menzinger, 1999)
2.4
The matrix determines the evolution of small perturbations from the limit cycle along tangent and normal directions, respectively. Based on the definitions of the unitary matrix U given by equation 2.1 and the Jacobian matrix given by equation 2.3, it can be shown that
2.5
where (see appendix  A for detailed derivations)
2.6
where . The notation means that the perturbation applied along the direction produces its effect along the direction . For example, the coefficient represents the rate of convergence of a tangent perturbation measured along the tangent direction of the limit cycle, that is, the amount of phase resetting induced by a sudden tangential perturbation of a figurative point. In a type 1 neural oscillator, dominates the other stability coefficients, which significantly simplifies the PRC prediction by neglecting the contribution of normal relaxation to the temporal phase resetting. Such an approximation is equivalent to the strongly attractive limit cycle assumption of Ermentrout and Kopell (1984, 1991a, 1991b), Hansel, Mato, and Meunier (1995), and Hoppensteadt and Izhikevich (1997). The coefficient is the rate of convergence of a normal perturbation measured along the normal direction to the limit cycle. In all planar dynamical systems, is zero, which means that a perturbation tangent to the limit cycle never induces a normal displacement of the figurative point. In other words, tangent perturbations can only advance or delay the phase of the figurative point without changing the amplitude of oscillations. However, , which means that normal perturbations also induce tangent effects, that is, phase resetting. Although other studies (e.g., Ali & Menzinger, 1999) alluded to the possible relationship between the stability coefficients given by equation 2.5 and the PRC, to our knowledge, the geometry of the temporal (phase) resetting has never been explicitly given in terms of the limit cycle stability coefficients and the geometric-to-temporal mapping.

In order to find the local tangent, xt, and normal, xn, geometric displacements along the limit cycle , which determines the PRC, we used the following algorithm:

Step 1. The limit cycle solution must be obtained by either (1) numerical integration of equation 1.2, as in our first detailed example that uses Bonhoeffer–van der Pol (BvP) model, or (2) analytically whenever the dynamical system is integrable, as in our second detailed example that uses a Poincaré model.

Step 2. Once the unperturbed limit cycle solution is known, equation 2.6 gives the coefficients along .

Step 3. The initial conditions for the linear differential equations given by equation 2.4 must be estimated based on the actual shape of the perturbation. Traditionally the PRC measures the amount of phase advance or delay when a brief current pulse affects only the transmembrane voltage, that is, the perturbation affects only the first equation of the dynamical system, equation 1.2. In the (traditional) PRC approach, the actual perturbation takes the form with respect to the fixed reference frame, (Oprisan & Canavier, 2002). In planar dynamical systems, the tangent and normal projections of the perturbation, the initial conditions for equation 2.4, can be found using the unitary matrix given by equation 2.1:
2.7
where represents the onset of the perturbation.
Step 4. A partial geometric resetting is determined by solving the linear and homogeneous differential equation for normal perturbations. For any planar dynamical systems, and, therefore, , whose solution is
2.8
where, with the initial condition .
Step 5. The next step is to determine the geometric displacement along the tangent direction by integrating the corresponding linear differential equations for dxt/dt. To this end, let us assume for a moment that the normal-to-tangent coupling, , vanishes. In this case, the only contribution to the phase resetting would be determined by pure tangent perturbations governed by the linear and homogeneous differential equation , which has the solution
2.9
where with the initial condition . However, the differential equation for tangent perturbations, which determines the total amount of geometric phase resetting, contains both pure tangent perturbations through the coefficient given by equation 2.9 and normal perturbations through :
2.10
Using the solution of equation 2.8 for normal perturbations and the solution of equation 2.9 for pure tangent perturbations, we determined the temporal evolution of cumulative tangent perturbations by integrating equation 2.10,
2.11
where the contributions of the normal perturbations to the PRC are determined by .

Step 6. Finally, the total tangential geometric resetting given by equation 2.11 can be converted into a temporal resetting by dividing it for the local vector flow F.

## 3.  Results

### 3.1.  Geometric-to-Temporal Phase Resetting Map for a Nonintegrable Neural Model.

This step-by-step geometrical method for computing the PRC by mapping a geometric displacement into a temporal phase resetting could be applied to any limit cycle oscillator. As a first concrete application to neurobiology, we explicitly show detailed calculations for the two-dimensional FitzHugh reduction of the Hodgkin-Huxley equations, Bonhoeffer–van der Pol equations (BvP) (FitzHugh, 1961), in dimensionless form:
3.1
Linear stability analysis of equations 3.1 shows that the steady state (x1ss=0, x2ss=0) is unstable and leads to a stable limit cycle (see FitzHugh, 1961; Lecar & Nossal, 1971). The steady state is at the intersection of the two nullclines f1=0 (voltage nullcline, dashed line in Figure 3C) and f2=0 (slow activation, dashed-dotted line in Figure 3C).
Figure 3:

Local stability of the limit cycle to small perturbations. Tangent perturbations of the figurative points can induce significant tangent effects, that is, phase resetting, if the limit cycle is unstable to tangent perturbations () (A). Additional phase resetting comes from normal perturbations to the limit cycle that are poised to induce a tangent effect through the coupling coefficient (B). For all two-dimensional limit cycles, and (D) has no direct contribution to the phase resetting. The regions of limit cycle instability due to tangent perturbations () coincide with the initial phases of the jumps between the depolarized and the hyperpolarized branches of the limit cycle. Normal perturbations induce phase resetting mostly along the depolarized and the hyperpolarized branches of the limit cycle () with additional narrow perturbation-sensitive regions at the crossing of the slow nullcline f2=0. The solid rhombuses mark the values of the stability coefficients at phase .

Figure 3:

Local stability of the limit cycle to small perturbations. Tangent perturbations of the figurative points can induce significant tangent effects, that is, phase resetting, if the limit cycle is unstable to tangent perturbations () (A). Additional phase resetting comes from normal perturbations to the limit cycle that are poised to induce a tangent effect through the coupling coefficient (B). For all two-dimensional limit cycles, and (D) has no direct contribution to the phase resetting. The regions of limit cycle instability due to tangent perturbations () coincide with the initial phases of the jumps between the depolarized and the hyperpolarized branches of the limit cycle. Normal perturbations induce phase resetting mostly along the depolarized and the hyperpolarized branches of the limit cycle () with additional narrow perturbation-sensitive regions at the crossing of the slow nullcline f2=0. The solid rhombuses mark the values of the stability coefficients at phase .

Based on the numerically obtained limit cycle solution (step 1 in section 2) and the unitary transform U from the fixed x=(x1, x2) to the moving reference frame given by equation 2.5, we computed the matrix B, which determines the stability of a limit cycle oscillator with respect to the normal and the tangent perturbations at different phases (step 2 in section 2). A negative value of any coefficient means that the perturbation decays exponentially along the respective direction; that is, the limit cycle is stable to perturbations applied at that phase. However, there are some phases along the limit cycle where perturbations have significant effects; the limit cycle is unstable at those phases. Since we are concerned with the relationship between the limit cycle stability to local perturbations measured by the coefficients and the PRC, our goal is to estimate the effect of tangent perturbations (to the unperturbed limit cycle) and convert (map) the geometric resetting into advances or delays of the temporal phase . We accounted for both the effect of , an obvious candidate to phase resetting contribution, and the coupling between a normal perturbation that can induce tangent effects because of the coupling coefficient in equation 2.5.

For the BvP model given by equations 3.1, the unstable regions of the limit cycle are shown in Figure 3. Tangent perturbations produce a significant phase resetting if (see Figure 3A), which happens only during the initial phases of the smooth transition from the depolarized to the hyperpolarized branch or vice versa (see the red thick lines in Figure 3C). The limit cycle is also sensitive to normal perturbations (see Figure 3D), during almost the same phases as for the above-discussed tangent perturbations (see the green thick lines in Figure 3C). However, does not induce any geometric phase resetting by itself; it needs a nonzero coupling term between the normal and the tangent modes in order to produce any phase resetting (see Figure 3B). The overlap of all tangent and normal instabilities leads to extreme sensitivity to perturbations and corresponds to the maxima of the PRC shown in Figure 6. In Figure 3 and all subsequent figures of this section, we marked with a solid rhombus a particular phase () that we followed throughout the entire section in order to explain how the PRC shown in Figure 6 relates to the stability coefficients shown in Figure 3. As seen from Figure 3, at all coefficients are negative: , , and .

Depending on the sign of the stability coefficient (see Figure 3C), the initial perturbation xn0 (the dashed line in Figure 4A) can exponentially decrease or increase during the time interval (see the solid line in Figure 4A). The elapsed time over which we followed the perturbations was between Pi/100 and Pi. Figure 4 and all subsequent figures in this section show representative examples with . For example, the amplitude of the normal perturbation xn at phase (see Figure 4A) is proportional to the initial normal projection of the perturbation, (see the empty rhombus in Figure 4A), multiplied by the exponential of the average during the interval (see equation 2.8). For , the average value of (step 4 in section 1). As a result, (see the solid rhombus in Figure 5A), which leads to a normal perturbation (see the solid rhombus in Figure 4A).

Figure 4:

Temporal evolution of perturbations. The evolution of normal perturbations (A) depends on the sign of . An initial normal perturbation xn0 (dashed line) decreases (if ) or increases (if ) over a time interval (solid line in panel A). Similarly, panel B shows the temporal evolution of a tangent perturbation xt0 (dashed line) that after shrinks (if ) and expands where (solid line). The open rhombus marks the value of the initial perturbation at phase , and the solid rhombus indicates the magnitude of the perturbation at with . (C) The total geometrical displacement includes the tangent displacement from panel B plus a scaled contribution of the normal perturbation from panel A that is determined by the cross-over coefficient (see equation 2.11).

Figure 4:

Temporal evolution of perturbations. The evolution of normal perturbations (A) depends on the sign of . An initial normal perturbation xn0 (dashed line) decreases (if ) or increases (if ) over a time interval (solid line in panel A). Similarly, panel B shows the temporal evolution of a tangent perturbation xt0 (dashed line) that after shrinks (if ) and expands where (solid line). The open rhombus marks the value of the initial perturbation at phase , and the solid rhombus indicates the magnitude of the perturbation at with . (C) The total geometrical displacement includes the tangent displacement from panel B plus a scaled contribution of the normal perturbation from panel A that is determined by the cross-over coefficient (see equation 2.11).

According to the step 5 in section 2, the amplitude of the tangent perturbation xt at phase (see the solid line in Figure 4B) is proportional to the initial tangent projection of the perturbation, (see the empty rhombus in Figure 4B), multiplied by the exponential of the average during the interval (see equation 2.9). For , the average value of . Notice that in agreement with Figure 3A, although , its average over 0.4Pi is positive because the average included both the small negative lobe and the large positive lobe between and . As a result, (see the solid rhombus in Figure 5B), which leads to a tangent perturbation (see the solid rhombus in Figure 4B).

Figure 5:

Temporal evolution of unit perturbations determined by (A), (B), and (C). Normal perturbations that induce normal displacements along the limit cycle (A) have a neutral region given by , which corresponds to . As noticed from panel A, the limit cycle is stable to normal perturbations since for . Similarly, tangent perturbations that induce tangent effects (phase resetting) also have neutral regions where , which lead to (B). For the cross-term , the neutral stability condition is (C).

Figure 5:

Temporal evolution of unit perturbations determined by (A), (B), and (C). Normal perturbations that induce normal displacements along the limit cycle (A) have a neutral region given by , which corresponds to . As noticed from panel A, the limit cycle is stable to normal perturbations since for . Similarly, tangent perturbations that induce tangent effects (phase resetting) also have neutral regions where , which lead to (B). For the cross-term , the neutral stability condition is (C).

Although the local divergences (instabilities) are exponential, when they are averaged over a wide range of phases, their effect is no longer dramatic. As expected from Figure 3A, the integral given by equation 2.9, , can be larger than unit and leads to local growth of perturbations (see Figure 5B), as we found in the case of a perturbation applied at (see Figure 4B). Since for most phases (see Figure 3D), the corresponding integral is always less than unit (see Figure 5A) and leads to an exponential decrease of any normal perturbation (see Figure 4A). The contribution of the cross-term is given by the integral Itn(t0, t) (see Figure 5C), which has both positive and negative values. We already computed the contribution of the pure tangent perturbations to the total geometric displacement and found that xt=−0.0014 (see the solid rhombus in Figure 4B). To this value, we need to add the contribution of the crossover term that takes a normal displacement (see the empty rhombus in Figure 4A) and converts it into a tangent displacement (i.e., multiply by the value 0.322 of the integral shown in Figure 5C at phase ). As a result, we find an additional tangent displacement of , which together with the previously computed pure tangent resetting gives us a total geometric displacement tangent to the limit cycle of −0.0046 (see the solid rhombus in Figure 4C).

Finally, based on step 6 in section 2, the total tangent geometric displacement given by equation 2.11 and shown in Figure 4C was converted to a temporal (or phase) advance/delay by dividing it to the local velocity of the flow given by equation 1.2.

Our predicted type 2 PRC based entirely on the geometrical reasoning is shown in Figure 6 with a dotted line. This PRC represents the phase shift induced at different phases around the limit cycle by a brief injected current . Some authors call this type of PRC an infinitesimal PRC (iPRC) to emphasize that this is the response of a neurobiological model to an arbitrarily brief current perturbation (Ermentrout, 1996; Gutkin, Ermentrout, & Reyes, 2005).

Figure 6:

Geometric-based prediction method for the PRC. The PRC obtained by numerically integrating the differential equations of BvP model given by equations 3.1 subject to a current pulse of duration (solid line) shows a typical type 2 resetting that includes both advances and delays of the subsequent spike. The geometrically based predicted PRC uses the numerically generated limit cycle solution of equations 3.1 and the stability coefficients given by equation 2.6 (dotted line). The improvement obtained by employing a moving reference frame and the coefficients is significant compared to previous approaches to PRC prediction for type 2 excitability class (see Figure 1C).

Figure 6:

Geometric-based prediction method for the PRC. The PRC obtained by numerically integrating the differential equations of BvP model given by equations 3.1 subject to a current pulse of duration (solid line) shows a typical type 2 resetting that includes both advances and delays of the subsequent spike. The geometrically based predicted PRC uses the numerically generated limit cycle solution of equations 3.1 and the stability coefficients given by equation 2.6 (dotted line). The improvement obtained by employing a moving reference frame and the coefficients is significant compared to previous approaches to PRC prediction for type 2 excitability class (see Figure 1C).

In any computer simulations, the current perturbation, however brief, is always finite. As a result, inevitably any numerically or experimentally generated PRC that was obtained by applying a finite duration perturbation to a neurobiological oscillator is not an iPRC. Numerically it is possible to shorten the duration of the perturbation in order to infer the shape of an iPRC or use the analytical method of the adjoint operator (see Ermentrout & Kopell, 1984, 1985, 1991a, 1991b; Ermentrout, 2002, for extensive examples). The comparison between the numerically generated PRC (obtained by perturbing the membrane voltage equations with a brief—0.5% of intrinsic period—rectangular current pulse) and the predicted iPRC based on the stability coefficients and equation shows a good match (Figure 6).

### 3.2.  Geometric-to-Temporal Phase Resetting Maps for an Integrable Neural Model.

Recently Ermentrout, Glass, and Oldeman (2012) used a Poincaré oscillator model with an SNIC bifurcation that has an analytic limit cycle solution that produces both unimodal type 1 and bimodal type 2 PRCs. We applied the geometric method to gain further qualitative insight into the relationship between the geometry of the limit cycle and the PRC. The model's equations are:
3.2
where parameter a indicates how strongly the limit cycle attracts nearby phase-space trajectories, b controls the distance to the SNIC bifurcation, and determines where the dynamics slow down (see Ermentrout et al., 2012). The model given by equations 3.2 has a SNIC bifurcation when b = 1 and a stable attracting limit cycle oscillations for . As Ermentrout et al. (2012) showed, this model has a limit cycle given by where and a period of oscillations (see step 1 in section 2).
Based on the definitions of the stability exponents given by equations 2.6 and 3.2, we found that (see step 2 in section 2)
3.3
As expected, only in equation 3.3 is affected by how strongly the limit cycle attracts nearby phase-space trajectories, and the other two coefficients that directly contribute to phase resetting depend on only the distance from the SNIC bifurcation through the parameter b (see Figure 7A). The geometric displacement from the unperturbed limit cycle is determined both by the coefficients given by equation 3.3 and the projection of the perturbation along the tangent and the normal directions to the limit cycle (see Figure 7B) given by (see step 3 in section 2)
3.4
The total contribution of the normal perturbations to the PRC (see step 4 in section 2), xn(t)=xn0Itt(t0, t)Itn(t0, t), was numerically computed (see Figure 7D). The contribution of tangent perturbations to the PRC (see step 5 in section 2) is determined by xt(t)=xt0Itt(t0, t), where , which is shown in Figure 7C for the same values of the parameter as in Ermentrout et al. (2012). We also found a closed form of the integral that accounts for the contribution of , , although an analytical expression for Itn(t0, t) could not be obtained.
Figure 7:

Normal and tangent perturbations for an integrable neural oscillator. (A) The tangent and normal stability coefficients of the integrable system given by equations 3.2 for b=0.9. The symmetrical shape of with both positive and negative values suggests that for model parameters that lead to significant normal displacements to the limit cycle, the resultant PRC is bimodal. (B) The tangent xt0 (solid line) and normal xn0 (dashed line) projections of the unit perturbation along the limit cycle for b=0.9 and . The tangent xt (C) and normal xn (D) contribution of a unit perturbation to the PRC for b=0.99 and a=10.0 and the same values of as in Ermentrout et al. (2012). The tangent contribution is significantly larger than the normal one for the above parameters.

Figure 7:

Normal and tangent perturbations for an integrable neural oscillator. (A) The tangent and normal stability coefficients of the integrable system given by equations 3.2 for b=0.9. The symmetrical shape of with both positive and negative values suggests that for model parameters that lead to significant normal displacements to the limit cycle, the resultant PRC is bimodal. (B) The tangent xt0 (solid line) and normal xn0 (dashed line) projections of the unit perturbation along the limit cycle for b=0.9 and . The tangent xt (C) and normal xn (D) contribution of a unit perturbation to the PRC for b=0.99 and a=10.0 and the same values of as in Ermentrout et al. (2012). The tangent contribution is significantly larger than the normal one for the above parameters.

As we notice from Figures 7C and 7D, for this particular dynamical system, the contribution of the normal resetting to the PRC through and is about 20 times smaller than that of the tangent perturbations determined by .

Finally, we compared the PRCs obtained in Ermentrout et al. (2012) by integrating equations 2.11 (see the traces marked with different geometric symbols in Figure 8) against our geometric-based predicted PRC (see the corresponding continuous lines in Figure 8) and found a very good match.

Figure 8:

The PRC of an integrable dynamical system. (A) The PRCs obtained by numerical integration of equations 3.2 with b=0.99, a=10.0, and (filled triangles), (filled circles), (filled squares), (filled rhombuses), and (filled inverted triangles) (adapted from Ermentrout et al., 2012). Geometric-based predictions of the PRCs (see section 2) are shown with solid lines and they overlap with data from Ermentrout et al. (2012).

Figure 8:

The PRC of an integrable dynamical system. (A) The PRCs obtained by numerical integration of equations 3.2 with b=0.99, a=10.0, and (filled triangles), (filled circles), (filled squares), (filled rhombuses), and (filled inverted triangles) (adapted from Ermentrout et al., 2012). Geometric-based predictions of the PRCs (see section 2) are shown with solid lines and they overlap with data from Ermentrout et al. (2012).

## 4.  Conclusion

In this letter, we have used an unambiguous definition of the geometrical (distance-based) phase that applies to every nonlinear oscillator even in high-dimensional spaces. We showed that external perturbations of limit cycle oscillators lead to local changes in the phase-space trajectory that can be subsequently converted into temporal advances or delays of the subsequent spikes. In order to measure the geometric resetting, we used a moving reference frame attached to the figurative point. The local eigenvalues on the limit cycle determined the relaxation speeds in the directions normal and tangent to the limit cycle at the point of the perturbation.

We previously found good agreement between the geometrically predicted PRC and the one obtained by direct integration of the model's equations only for the type 1 excitability class (Oprisan & Canavier, 2002, 2006). Therefore, our main concern in this letter was to improve the prediction power of the geometrical method for the type 2 excitability class.

For this purpose, we used both a nonintegrable type 2 neural oscillator model (Bonhoeffer-van der Pol) and an integrable (Poincaré) model in order to compare the predicted PRC based on the geometric-to-temporal mapping against the numerically generated PRCs. We found good agreement between the geometrically based PRC prediction and the ones obtained by integrating the corresponding differential equation subject to brief current perturbations.

The improvement in the prediction accuracy of type 2 PRC is obvious when comparing Figure 1C against Figure 6 or Figure 8. In our opinion, the main source of error in applying this geometrical method that uses the stability coefficients to trace the tangent (phase resetting) effect of normal perturbations specific to a type 2 excitable membrane is due to the linearization of the perturbation equation. Perhaps a better approximation of the perturbation equation that would include quadratic terms could reduce the estimation error. However, a drawback of such an incremental approach is that the geometric method that we presented here would lose its immediate and intuitive, geometrical interpretation.

Another advantage offered by this geometrical approach to the PRC prediction is that it can accommodate arbitrarily shaped stimulus current simply by substituting the particular time-dependent stimulus in equation 2.7. As Guillamon and Huguet (2009) noted, in the (traditional) theory of phase resetting, we are interested in only the effect of the perturbations along the direction of the membrane voltage, which gives the PRC. This geometrical method can be extended to accommodate multidimensional perturbations.

Can the approach presented here be expanded to higher dimensions? We tested the method with the three-dimensional Hindmarsh- Rose model (Hindmarsh & Rose, 1984) with the tangent vector (as we suggested in section 2) and different orthogonal variations of the normal vectors. Regardless of the particular selection of normal vectors, we found that the stability exponent always vanishes. Furthermore, by appropriately selecting the two normal vectors in the unitary matrix, it is possible to cancel out at least one more of the two nonzero subdiagonal elements of the matrix B. We are currently searching for a change of variables that would allow us to write a simple closed form for the PRC similar to the two-dimensional case.

Can this method be extended to physiological neurons? Although we do not have access a priori to the equations that govern the electrical activity of a given neuron, this method can still apply to physiological neurons because we have observed (Oprisan et al., 2003) that the geometrical-to-temporal phase relationship that we obtained using the model's equations is similar to that obtained from membrane potential records after delay-embedding reconstruction of the phase's space attractor. More generally, we can replace the described Euclidean-based approach (which requires analytic expressions of the vector field) by a computational geometrical phase construction based entirely on a delayed embedding using the fast variable record (Takens, 1981). In order to obtain the slow variable(s), filtering techniques may be required (see a detailed phase-space reconstruction of pyloric attractor from experimental data in Oprisan et al., 2003). Once we have an appropriate phase-plane reconstruction, the geometrical-to-temporal phase mapping can be obtained and then the predicted phase resetting induced by external perturbations. We predict that these two steps can give us the PRC for neurons (such as cortical neurons) that are type 1 oscillators (Wilson & Bower, 1989; Traub & Miles, 1991). Experimentally obtained PRCs can then be used to confirm the validity of our theoretical analysis.

The geometrical method described here applies to perturbations that produce local displacements from the limit cycle and allow the figurative point to relax back to the unperturbed limit cycle as shown in Figure 2. Although the perturbations we considered are not necessarily weak or infinitesimal, this local geometrical method assumes that the phase resetting is entirely determined by the local eigenvalues on the limit cycle that dictate different relaxation speeds in the directions normal and tangent to the limit cycle at the point of the perturbation. In contrast, we showed that very strong inhibitory perturbations driving the pyloric dilator (PD) neuron in the stomatogastric ganglion of Homarus americanus force the figurative point to switch between the depolarized (bursting) branch of the limit cycle and the hyperpolarized (silent) one (see Oprisan et al., 2003 and the extension to excitatory perturbations by Maran et al., 2011). In order to correctly predict the PRC in response to such very strong perturbations that produce global changes in the location of the figurative point, an additional geometrical mapping is required (Oprisan et al., 2003). This global (or “normal” to the separatrix between the depolarized and hyperpolarized branches of the limit cycle) map was obtained from the reconstructed limit cycle based on membrane potential recordings by using the delay-embedding method (see Figure 3A in Oprisan et al., 2003). The global geometrical phase mapping (see Figure 3B in Oprisan et al., 2003) paired phases on the depolarized and hyperpolarized branches of the limit cycle. “The normal component of phase resetting is given by on the depolarized branch, and by on the hyperpolarized branch.” When the limit cycle perturbations are very strong (Oprisan et al., 2003; Maran et al., 2011), such a global phase-resetting approach (or “normal component” to the separatrix between the depolarized and hyperpolarized branches of the limit cycle) is required, in addition to the local analysis of phase-space relaxation detailed in this letter.

Finally, there are similarities between our local method of predicting phase resetting and those given in the appendix of Ermentrout and Kopell (1991a, 1991b), but there are some significant differences as well (see Oprisan & Canavier (2002, 2006) and appendix  B in this letter for detailed comparisons with similar mathematical approaches). Although previous work contains expressions that are mathematically similar to the ones presented in this letter, they do not rely on the geometrical intuition presented here. This intuition allows us to extend the PRC analysis to type 2 model neurons and even to physiological neurons for which we do not know the equations of state but for which we can reconstruct the phase-space attractor.

## Appendix A: Stability Coefficients of a Planar Limit Cycle in a Moving Reference Frame

For small perturbations, , from the unperturbed limit cycle , equation 1.2 reduces to
A.1
The unitary matrix U that transforms the coordinates of the fixed reference frame x=(x1, x2)T to the moving reference frame , is
A.2
As a result, equation A.1, which describes the evolution of limit cycle perturbations in x=(x1, x2)T, rewrites in coordinates as follows:
A.3
where
A.4
As a result, is given by:
A.5
and UAU−1(f21+f22)2 is given by
A.6
Based on equations A.5 and A.6, the matrix that determines the evolution of small perturbations from the limit cycle along both the tangent and that normal directions is
A.7

## Appendix B: Comparison with Previous Studies

There are similarities between our method of predicting phase resetting and those given in the appendix of Ermentrout and Kopell (1991a, 1991b), but there are some significant differences as well. Our derivation uses geometrical intuition, whereas that of Ermentrout and Kopell is analytical. In Ermentrout and Kopell (1991a, 1991b), a system of N state variables describes neural oscillators using equations similar to equation 1.2. The original N variables were changed to a more convenient phase variable , which parameterizes the limit cycle solution of an isolated oscillator, plus N−1 variables that describe the normal displacements from the limit cycle. Their approach is very similar to our geometric approach of changing the variables to the tangent and the normal coordinates with respect to the limit cycle. However, in Ermentrout and Kopell (1991a, 1991b), the N−1 variables can be ignored by assuming a strong (“infinitely”), attractive limit cycle that renders these N−1 variables negligibly small. Furthermore, Ermentrout and Kopell assumed that the phase variables of (two) coupled oscillators are never perturbed independently, which further reduced the complexity of the dynamical system to a single-phase variable.

Both our derivation and that of Ermentrout and Kopell take advantage of the fact that the perturbation is in one variable (membrane voltage) only. For example, they assumed the following coupling function in a system of two coupled Morris-Lecar oscillators , which is equivalent to the scaled instantaneous current perturbation from step 3 in section 2 of this letter. In Ermentrout and Kopell, the evolution equation for the reduced angular phase is , where i is the instantaneous (synaptic) current and is the constant angular frequency of the limit cycle oscillator. A similar expression for appears in Hoppensteadt and Izhikevich (1997). As explained in Ermentrout and Kopell (1984), in case is not constant, for example, , it is always possible to change variables, for example, where with . The second term in their formulate of , , simply gives the phase resetting due to a perturbation i by projecting the voltage perturbation along the tangent direction to the limit cycle. While our geometrical approach has the same goal of estimating the PRC, the geometrical intuition behind this reasoning is described only in this letter.

Assuming for a moment that only pure tangent displacements to the limit cycle are considered when predicting the PRC, we showed in this paper that (see equation 2.9) xt(t)=xt0Ttt(t0, t) where . During step 6 of in section 2, this geometrical tangent perturbation of the limit cycle due to a voltage perturbation was converted to a temporal phase resetting by dividing xt(t) to the local vector flow amplitude. Therefore, in our geometrical approach, the phase resetting due solely to a tangent perturbation is proportional to . We notice that is the adjoint, or the infinitesimal PRC, as given by equation 3.2 of Hansel et al. (1995) and by equation A.4 of Ermentrout and Kopell (1991b). We must emphasize, however, that this immediate correspondence between our geometrical approach and those in Ermentrout and Kopell (1984, 1991a, 1991b), Hansel et al. (1995), and Hoppensteadt and Izhikevich (1997) holds providing that perturbations normal to the limit cycle are ignored. Although previous work in the literature contains expressions that are mathematically similar to the ones presented in this letter, they do not contain the geometrical intuition presented here.

## Acknowledgments

This work was supported by the National Science Foundation IOS CAREER award 1054914 to me. I am grateful to the editor and to the two anonymous referees for their generous and insightful comments that led to a substantial improvement of the letter.

## References

Abramovich-Sivan
,
S.
, &
Akselrod
,
S.
(
1998a
).
A phase response curve based model: Effect of vagal and sympathetic stimulation and interaction on a pacemaker cell
.
Journal of Theoretical BioIogy
,
192
,
567
579
.
Abramovich-Sivan
,
S.
, &
Akselrod
,
S.
(
1998b
).
A PRC based model of a pacemaker cell: Effect of vagal activity and investigation of the respiratory sinus arrhythmia
.
Journal of Theoretical BioIogy
,
192
,
219
243
.
Ali
,
F.
, &
Menzinger
,
M.
(
1999
).
On the local stability of limit cycles
.
Chaos
,
9
(
2
),
348
356
.
Alving
,
B. O.
(
1968
).
Spontaneous activity in isolated somata of Aplysia pacemaker neurons
.
Journal of General Physiology
,
51
,
29
45
.
Baxter
,
R.
,
Canavier
,
C.
,
Lechner
,
H.
,
Butera
,
R.
,
Clark
,
A.D.J.
, &
Byrne
,
J.
(
2000
).
Coexisting stable oscillatory states in single cell and multicellular neuronal oscillators
. In D. Levine, V. Brown, & T. Shirley (Eds.),
Oscillations in neural systems
.
Hillside, NJ
:
Erlbaum
.
Bellman
,
R.
(
1997
).
Introduction to matrix analysis
.
:
Society for Industrial and Applied Mathematics
.
Brown
,
E.
,
Moehlis
,
J.
, &
Holmes
,
P.
(
2004
).
On the phase reduction and response dynamics of neural oscillator populations
.
Neural Comput.
,
16
(
4
),
673
715
.
Cheron
,
G.
,
Cebolla
,
A.
,
Saedeleer
,
C.
,
Bengoetxea
,
A.
,
Leurs
,
F.
,
Leroy
,
A.
, et al
(
2007
).
Pure phaselocking of beta/gamma oscillation contributes to the n30 frontal component of somatosensory evoked potentials
.
BMC Neuroscience
,
8
(
75
),
1
11
.
Dror
,
R.
,
Canavier
,
C.
,
Butera
,
R.
, &
Byrne
,
J.C.J.
(
1999
).
A mathematical criterion based on phase response curves for the stability of a ring network of oscillators
.
Biological Cybernetics
,
80
,
11
23
.
Ermentrout
,
G. B.
(
1994
).
Neural modelling and neural networks
.
Oxford
:
Pergamon Press
.
Ermentrout
,
G. B.
(
1996
).
Type l membranes, phase resetting curves, and synchrony
.
Neural Computation
,
8
(
5
),
979
1001
.
Ermentrout
,
G. B.
(
2002
).
Simulating, analyzing, and animating dynamical systems: A guide to XPPAUT for researchers and students, 14
.
:
Society for Industrial and Applied Mathematics
.
Ermentrout
,
G. B.
,
Glass
,
L.
, &
Oldeman
,
B. E.
(
2012
).
The shape of phase-resetting curves in oscillators with a saddle node on an invariant circle bifurcation
.
Neural Computation
,
24
,
3111
3125
.
Ermentrout
,
G. B.
, &
Kopell
,
N.
(
1984
).
Frequency plateaus in a chain of weakly coupled oscillators, I
.
SIAM Journal on Mathematical Analysis
,
15
,
215
237
.
Ermentrout
,
G. B.
, &
Kopell
,
N.
(
1985
).
The behavior of rings of coupled oscillators
.
Journal of Mathematical Biology
,
23
,
55
74
.
Ermentrout
,
G. B.
, &
Kopell
,
N.
(
1986
).
Parabolic bursting in an excitable system coupled with a slow oscillation
.
SIAM Journal on Applied Mathematics
,
46
,
233
253
.
Ermentrout
,
G. B.
, &
Kopell
,
N.
(
1991a
).
Multiple pulse interactions and averaging in systems of coupled neural oscillators
.
Journal of Mathematical Biology
,
29
,
195
217
.
Ermentrout
,
G. B.
, &
Kopell
,
N.
(
1991b
).
Oscillator death in systems of coupled neural oscillators
.
SIAM Journal on Applied Mathematics
,
50
,
125
146
.
FitzHugh
,
R.
(
1961
).
Impulses and physiological states in theoretical models of nerve membrane
.
Biophysical Journal
,
1
,
445
466
.
Guillamon
,
A.
, &
Huguet
,
G.
(
2009
).
A computational and geometric approach to phase resetting curves and surfaces
.
SIAM Journal on Applied Dynamical Systems
,
8
(
3
),
1005
1042
.
Gutkin
,
B. S.
,
Ermentrout
,
G. B.
, &
Reyes
,
A. D.
(
2005
).
Phase-response curves give the responses of neurons to transient inputs
.
Journal of Neurophysiology
,
94
(
2
),
1623
1635
.
Haenschel
,
C.
,
Linden
,
D.
,
Bittner
,
R.
,
Singer
,
W.
, &
Hanslmayr
,
S.
(
2010
).
Alpha phase locking predicts residual working memory performance in schizophrenia
.
Biological Psychiatry
,
68
,
595
598
.
Hansel
,
D.
,
Mato
,
G.
, &
Meunier
,
C.
(
1995
).
Synchronization in excitatory neural networks
.
Neural Computation
,
7
,
307
337
.
Hindmarsh
,
J.
, &
Rose
,
R.
(
1984
).
A model of neuronal bursting using three coupled first order differential equations
.
Proceedings of the Royal Society
,
B 221
,
87
102
.
,
F. C.
, &
Izhikevich
,
E. M.
(
1997
).
Weakly connected neural networks
.
New York
:
Springer-Verlag
.
Izhikevich
,
E. M.
(
2000
).
Phase equations for relaxation oscillators
.
SIAM Journal on Applied Mathematics
,
60
,
1789
1804
.
Izhikevich
,
E. M.
(
2007
).
Dynamical systems in neuroscience: The geometry of excitability and bursting
.
Cambridge, MA
:
MIT Press
.
Lecar
,
H.
, &
Nossal
,
R.
(
1971
).
Theory of threshold fluctuations in nerves II. Analysis of various sources of membrane noise
.
Biophysical Journal
,
11
(
12
),
1068
1084
.
Maran
,
S. K.
,
Sieling
,
F. H.
,
Demla
,
K.
,
Prinz
,
A. A.
, &
Canavier
,
C. C.
(
2011
).
Responses of a bursting pacemaker to excitation reveal spatial segregation between bursting and spiking mechanisms
.
Journal of Computational Neuroscience
,
31
,
419
440
.
Murray
,
J. D.
(
1993
).
Mathematical biology
(2nd ed.).
New York
:
Springer
.
Oprisan
,
S.
(
2009
).
Stability of synchronous oscillations in a periodic network
.
International Journal of Neuroscience
,
119
,
482
491
.
Oprisan
,
S. A.
(
2012a
).
Existence and stability criteria for phase-locked Modes in ring networks using phase-resetting curves and spike time resetting curves
.
New York
:
Springer
.
Oprisan
,
S. A.
(
2012b
).
A geometric approach to phase resetting estimation based on mapping temporal to geometric phase
.
New York
:
Springer
.
Oprisan
,
S.
, &
Boutan
,
C.
(
2008
).
Prediction of entrainment and 1:1 phase-locked modes in two-neuron networks based on the phase resetting curve method
.
International Journal of Neuroscience
,
118
,
867
890
.
Oprisan
,
S. A.
, &
Canavier
,
C. C.
(
2002
).
The influence of limit cycle topology on the phase resetting curve
.
Neural Computation
,
14
,
1027
1057
.
Oprisan
,
S. A.
, &
Canavier
,
C. C.
(
2006
).
The structure of instantaneous phase resetting in a neural oscillator
. In A. A. Minai & Y. Bar-Yam (Eds.),
Unifying themes in complex systems
(pp.
223–231
).
Berlin
:
Springer
.
Oprisan
,
S. A.
,
Prinz
,
A. A.
, &
Canavier
,
C. C.
(
2004
).
Phase resetting and phase locking in hybrid circuits of one model and one biological neuron
.
Biophysical Journal
,
87
(
4
),
2283
2298
.
Oprisan
,
S. A.
,
Thirumalai
,
V.
, &
Canavier
,
C. C.
(
2003
).
Dynamics from a time series: Can we extract the phase resetting curve from a time series
?
Biophysical Journal
,
84
(
5
),
2919
2928
.
Pinsker
,
H.
(
1977
).
Aplysia bursting neurons as endogenous oscillators: I: Phase-response curves for pulsed inhibitory synaptic input
.
Journal of Neurophysiology
,
40
,
527
543
.
Rutishauser
,
U.
,
Ross
,
I.
,
Mamelak
,
A.
, &
Schuman
,
E.
(
2010
).
Human memory strength is predicted by theta-frequency phase-locking of single neurons
.
Nature
,
464
,
903
907
.
Shores
,
T. S.
(
2007
).
Applied linear algebra and matrix analysis
.
New York
:
Springer
.
Takens
,
F.
(
1981
).
Detecting strange attractors in turbulence
. In D. Rand & L. Young (Eds.),
Proceedings of the Symposium on Dynamical Systems and Turbulence
.
Warwick
:
University of Warwick
.
Traub
,
R.
, &
Miles
,
R.
(
1991
).
Neural networks of the hippocampus
.
Cambridge
:
Cambridge University Press
.
Wilson
,
M.
, &
Bower
,
J.
(
1989
).
The simulation of large scale neural networks
.
Cambridge, MA
:
MIT Press
.
Winfree
,
A. T.
(
1980
).
The geometry of biological time
.
New York
:
Springer
.
Wingerden
,
M.
,
Vinck
,
M.
,
Lankelma
,
J.
, &
Pennartz
,
C.
(
2010
).
Learning-associated gamma-band phase-locking of action? Outcome selective neurons in orbitofrontal cortex
.
Journal of Neuroscience
,
30
,
10025
10038
.
Zacksenhouse
,
M.
, &
Ahissar
,
E.
(
2006
).
Temporal decoding by phase-locked loops: Unique features of circuit-level implementations and their significance for vibrissal information processing
.
Trends in Neural Computation
,
18
,
1611
1636
.