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

*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,

*P*: . We called the

_{i}*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 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.

**F**: 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,

*x*

_{1}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*+

*P*) that attracts all points in a neighborhood of as .

_{i}As a concrete example, we show the phase-plane (*x*_{1}, *x*_{2}) 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 *P _{i}* 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 *C _{m}* is the membrane capacitance, is the duration of the stimulus,

*f*

_{1}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 **v _{i}** 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.

**v**. Since we selected clockwise motion around the limit cycle (see Figure 2), the normal vector

_{2}**v**is given by . As a result, the unitary matrix is

_{2}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).

In order to find the local tangent, *x _{t}*, and normal,

*x*, geometric displacements along the limit cycle , which determines the PRC, we used the following algorithm:

_{n}**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: where represents the onset of the perturbation.

**Step 5.**The next step is to determine the geometric displacement along the tangent direction by integrating the corresponding linear differential equations for

*dx*/

_{t}*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 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 : 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, 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.

*x*

_{1ss}=0,

*x*

_{2ss}=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

*f*

_{1}=0 (voltage nullcline, dashed line in Figure 3C) and

*f*

_{2}=0 (slow activation, dashed-dotted line in Figure 3C).

Based on the numerically obtained limit cycle solution (step 1 in section 2) and the unitary transform *U* from the fixed **x**=(*x*_{1}, *x*_{2}) 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 *x*_{n0} (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 *P _{i}*/100 and

*P*. Figure 4 and all subsequent figures in this section show representative examples with . For example, the amplitude of the normal perturbation

_{i}*x*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).

_{n}According to the step 5 in section 2, the amplitude of the tangent perturbation *x _{t}* 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.4

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

_{i}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 *I _{tn}*(

*t*

_{0},

*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

*x*=−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).

_{t}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).

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.

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

*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) The total contribution of the normal perturbations to the PRC (see step 4 in section 2),

*x*(

_{n}*t*)=

*x*

_{n0}

*I*(

_{tt}*t*

_{0},

*t*)

*I*(

_{tn}*t*

_{0},

*t*), was numerically computed (see Figure 7D). The contribution of tangent perturbations to the PRC (see step 5 in section 2) is determined by

*x*(

_{t}*t*)=

*x*

_{t0}

*I*(

_{tt}*t*

_{0},

*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

*I*(

_{tn}*t*

_{0},

*t*) could not be obtained.

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 .

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

**x**=(

*x*

_{1},

*x*

_{2})

^{T}, rewrites in coordinates as follows: where

## 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) *x _{t}*(

*t*)=

*x*

_{t0}

*T*(

_{tt}*t*

_{0},

*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

*x*(

_{t}*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.