## Abstract

We introduce a simple two-dimensional model that extends the Poincaré oscillator so that the attracting limit cycle undergoes a saddle node bifurcation on an invariant circle (SNIC) for certain parameter values. Arbitrarily close to this bifurcation, the phase-resetting curve (PRC) continuously depends on parameters, where its shape can be not only primarily positive or primarily negative but also nearly sinusoidal. This example system shows that one must be careful inferring anything about the bifurcation structure of the oscillator from the shape of its PRC.

## 1. Introduction

Stimuli delivered to biological oscillators affect the timing of subsequent cycles. Because of the functional importance of this phase resetting, theoretical and experimental studies have probed resetting in diverse settings (Hoppensteadt & Keener, 1982; Winfree, 2000; Glass & Mackey, 1988; Izhikevich, 2007; Canavier, 2006; Rinzel & Ermentrout, 1998; Ermentrout & Terman, 2010). The results of such experiments are typically represented by phase-resetting curves (PRCs), which give the change of the phase of an oscillation as a function of the phase of a stimulus, and phase transition curves (PTCs), which give the new phase of an oscillation as a function of the phase of a stimulus (precise definitions are given below). If the original oscillation rapidly reestablishes itself following a pulsatile stimulus, then these curves can be used to analyze the effects of inputs to biological oscillations as well as the effects of coupling and noise on systems of biological oscillators. For example, phase-resetting curves can be used to predict the effects of periodic stimuli to oscillating heart cells (Guevara, Glass, & Shrier, 1981) and can also be used to analyze synchronization in systems of oscillators with pulsatile coupling (say, neurons with short-lasting synapses or fireflies) (Strogatz, 1994). The shape of the PRC determines the manner in which a pair of mutually connected oscillators will synchronize and how sensitive the locking will be to noise and heterogeneities (Hansel, Mato, & Meunier, 1995; Ermentrout, 1996). The shape also plays a crucial role in how oscillators react to correlated noise and other signals (Ermentrout, Galán, & Urban, 2008; Marella & Ermentrout, 2008). Thus, there is a good deal of interest in what aspects of an oscillator determine the shape of the PRC.

A common rule of thumb about PRCs in neuroscience is that oscillators that are near a saddle node on an invariant circle (SNIC) are of one sign (usually positive) (Ermentrout, 1996; Brown, Moehlis, & Holmes, 2004; Canavier, 2006), while those in which the oscillation arises from a Hopf bifurcation are bimodal (Ermentrout & Kopell, 1984; Brown et al., 2004). Recently Achuthan, Butera, and Canavier (2011) showed that as one moves away from the SNIC, a substantial bimodality developed in the PRC of the Morris-Lecar neural model. The main goal of this letter is to extend the results of a recent model (Krogh-Madsen, Butera, Ermentrout, & Glass, 2012), which demonstrates that PRCs that are sinusoidal in shape can also occur in systems that are arbitrarily close to a SNIC bifurcation. This analytically solvable system allows one to continuously transition from PRCs that are of one sign to PRCs that are sinusoidal.

## 2. Terminology

*T*and that we define a phase, , which is a time-like variable that increases at a constant rate of 1. At some phase , a brief perturbation is applied to one or more of the variables in the oscillator. This leads to a sequence of times, at which the oscillator crosses the zero phase line, which must be defined someplace. For excitable systems such as nerve cells or heart cells, it is common to take the point of maximal slope of the upstroke of an action potential as the point of zero phase, in which case it is called the fiducial point. Provided the state point following the stimulus is in the basin of attraction of the limit cycle, there will be an asymptotic return to the limit cycle in the limit . In the absence of the perturbation,

*T*=

_{j}*jT*, where

*j*is an integer. From the timing of the subsequent action potential upstrokes, we define the PTC mapping to the new phase, : For the situation in which the perturbation is small, we define the PRC, which is defined by asymptotic phase shift as In some situations that arise in experiments (e.g., transients following a stimulus in which several beats are skipped before an oscillation is reestablished), this definition of the PRC curve will not be applicable.

In this letter, we scale the phase so that , , and , but in other contexts, the phase here can be multiplied by or . The degree to which the value of *T _{j}*−

*T*

_{j−1},

*j*>1 differs from

*T*depends on the rate at which the rhythm returns to the limit cycle oscillation.

*isochron*(Guckenheimer, 1975; Winfree, 2000). Each point in the basin of attraction of a stable limit cycle will asymptotically approach the limit cycle as . Consider two trajectories in phase space

*x*(

_{i}*t*) and

*x*(

_{j}*t*) starting from points

*x*(0) and

_{i}*x*(0), where

_{j}*x*(0) is a point lying on the limit cycle and

_{i}*x*(0) is a point in the basin of attraction of the limit cycle. The points

_{j}*x*(0) and

_{i}*x*(0) are said to lie on the same isochron if the distance as . The isochrons foliate the basin of attraction of the limit cycle, so that each point in the basin of attraction of the limit cycle lies on an isochron.

_{j}## 3. A Simple Model of a Nonlinear Oscillator

There are several mechanisms by which a system of differential equations with a stable equilibrium point can lose the stability of this equilibrium and start to oscillate. Two of the best-known bifurcations are the supercritical Hopf bifurcation and the saddle node on an invariant circle (SNIC, also known as SNIPER for infinite-period) bifurcation. In the former case, the equilibrium exists but loses stability, and a small-amplitude limit cycle emerges. There appears to be widespread belief that for systems with a SNIC bifurcation, the PRC is nonnegative so that all stimuli will lead to an advance in the timing of subsequent action potentials, while near a Hopf bifurcation, the PRC will have a positive and negative lobe (Canavier, 2006). Our main point in this letter is to explain the shape of the PRC for a system with a SNIC bifurcation.

*a*and

*q*are positive. This equation (or similar equations with a circular limit cycle attractor) has been called the Poincaré oscillator (Glass & Mackey, 1988). In the limit when

*q*=0 or , the PRC and PTC can be readily computed using simple trigonometry, and the resulting equation has been called the radial isochron clock (Winfree, 2000; Hoppensteadt & Keener, 1982). This model has been useful for understanding the effects of periodic stimulation of the heart (Guevara & Glass, 1982).

*b*=1 and stable attracting limit cycle oscillations for . For

*b*>1, the limit cycle is lost, but there are two additional fixed points centered around . Thus, the parameter determines where the dynamics slows down. Figure 1 shows the geometry of this model. The limit cycle of this model is given by

*r*=1 and satisfying which can be solved by quadrature From these formulas, we find , , −

*T*/2<

*t*<

*T*/2, and defines the phase.

## 4. Phase Resetting

Understanding the shape of the PRCs in models is facilitated by understanding the geometry of the isochrons. As is well known (Glass & Mackey, 1988), for the Poincaré oscillator with *q*=0, the isochrons are equally spaced around the cycle and are aligned along the radii of the circle defined by the limit cycle oscillation. For the modified equation defined by equation 3.2 with 1>*b*, the isochrons accumulate in the neighborhood of the incipient saddle node bifurcation at .

*a*=1, 10,

*b*=0.99, 0.999, 0.9999. These isochrons were computed using AUTO (Doedel & Oldeman, 2012) via the algorithm described in Osinga and Moehlis (2010). This algorithm needs the Floquet eigenfunction corresponding to the nontrivial (in this case, stable) Floquet multiplier of the periodic orbit. Specifically, let

*dX*/

*dt*=

*F*(

*X*) be a differential equation in , and let

*U*(

*t*) be an exponentially stable limit cycle solution with period

*T*. Let

*A*(

*t*)≔

*D*(

_{X}F*U*(

*t*)) be the linearization of

*F*evaluated along the limit cycle. Then the Floquet eigenfunction

*V*(

*t*) can be defined as the unique

*T*-periodic solution to where

*V*(0) is the corresponding eigenvector for the time-

*T*map evaluated at

*t*=0. Moreover,

*V*(0) is the tangent vector to the isochron at that point. The algorithm collects starting points of time-

*nT*maps that, using boundary value methods, are forced to end at , where varies between 0 and and .

The Floquet multiplier becomes extremely small for large values of *a* and *b*, for example, for *a*=10, *b*=0.9999, which caused numerical difficulties when computed directly using the techniques described in Doedel, Kooi, van Voorn, and Kuznetsov (2008). We avoided these difficulties by first computing the Floquet eigenfunction for smaller parameter values: *a*=0.5 and *b*=0. Then the phase of the extended system consisting of the differential equations for the periodic orbit *U*(*t*) and its eigenfunction *V*(*t*) was shifted so that the desired isochron was at *t*=0. Finally, the extended system was continued in *a*, , and *T* and then in *b*, , and *T*, to reach the desired values of *a* and *b*. This homotopy approach, in combination with a high number of mesh points (the AUTO-constant NTST was set to 700), allowed us to reliably find the Floquet eigenfunction and then run the algorithm from Osinga and Moehlis (2010) that computes the actual isochrons.

From Figure 2, we note that as the value of *b* gets closer to 1, the isochrons are more tightly packed around , reflecting the slower angular velocity at that point, and as the value of *a* increases, the curvature of the isochrons decreases, as expected. For a different value of , the geometry can be found by an appropriate rotation of the images in Figure 2 so that the isochrons are always clustered in the region of .

As we have recently discussed in Krogh-Madsen et al. (2012), the shape of the PRC will depend on the location of the state point following a perturbation relative to and the fiducial point of the limit cycle oscillator. Thus, depending on the parameters, the PRC might be mostly positive, mostly negative, or bimodal. These situations are illustrated in Figure 3a, which shows the PRC for a stimulus applied to *u* with magnitude *s*=0.1 and several values of .

When , the PRC is almost entirely negative; when , the PRC is almost entirely positive; and when , the PRC is bimodal. The reason for this behavior is clear, and it results from the accumulation of the isochrons in the neighborhood of the incipient fixed point at . Thus, despite being very close to the saddle node bifurcation, the PRC can look very sinusoidal, and in particular, it is quite far from being strictly of one sign. To give the reader some sense of what the time traces of the modified oscillator look like, we plot *U*(*t*) for three values of in Figure 4 corresponding to the traces in Figure 3. The trace of is similar to the shape of an action potential.

The trajectories in Figure 4 suggest that the shape of the action potential might also be an indicator of whether the PRC has a negative lobe near the SNIC. As noted, for example, in Rinzel and Ermentrout (1998), near a SNIC, there is sometimes a small negative region after the peak of the action potential (compare the curve for in Figure 3). Geometrically, we would expect that if the slow part of the trajectory is at a place in the phase plane that is nearly opposite the peak of the potential, as, for example, is seen in the Hindmarsh-Rose model (Hindmarsh & Rose, 1984), then it may be possible to greatly increase the negative region of the PRC near the SNIC (compare the curve for in Figure 3).

### 4.1. Small-Amplitude Perturbations and Adjoints.

For many biological applications in which oscillators weakly interact, it is adequate to consider small-amplitude perturbations of the limit cycle. From a mathematical perspective, this limit is of interest because it is possible to carry out analytical study of the PRCs.

*s*is small, we expect that the PRC will be of the form , that is, it depends linearly on the amplitude

*s*. Here the function is the infinitesimal response to a stimulus applied to the

*n*th component of the limit cycle. Several authors (Izhikevich, 2007; Ermentrout & Terman, 2010) have shown that is a component in a certain vector quantity, the adjoint , which is closely related to the limit cycle. Indeed, is the gradient of the isochron map evaluated on the limit cycle (Izhikevich, 2007). The adjoint is easy to compute numerically for stable limit cycles. Specifically,

*Z*(

*t*) is the unique

*T*-periodic solution to where

*U*(

*t*) and

*A*(

*t*) are defined above at the definition of the Floquet eigenfunction in equation 4.1.

*Z*(

*t*) explicitly because of its close relationship with the original limit cycle. For example, consider the scalar differential equation on the circle, where

*f*(

*x*+1)=

*f*(

*x*) and

*f*(

*x*)>0 and continuous. Then there exists a solution

*U*(

*t*) defined implicitly by where

*u*(0)=0,

*u*(

*t*+

*T*)=

*u*(

*t*)+1 for all

*t*, and . The adjoint , and since

*f*(

*u*)>0, the adjoint is strictly positive. For with

*I*>0, the oscillator has a period and , which is strictly positive. Adjoints and PRCs that are nonnegative are called class I PRCs.

^{1}PRCs that have bimodal shapes (substantial regions where there is phase advance and phase delay) are called class II PRCs.

*a*and

*q*. Thus, for the simple Poincaré oscillator, the PRC is bimodal.

*r*=1 and are the solutions to the adjoint equation in polar form: We can readily solve for . Then we could solve the scalar equation for

*z*, where and

_{r}*h*(

*t*+

*T*)=

*h*(

*t*), and apply the appropriate normalization criteria. To begin, we use some simple approximations. We seek periodic solutions

*z*to equation 4.4. Let , and . Then we can rewrite equation 4.4 as

_{r}Now it follows from standard singular perturbation theory that (indeed, set , then *Y*(*t*)=−*h*(*t*) and , which is bounded as as it is independent of ). Thus, .

*u*component of the limit cycle as approximately If we choose , then with some trigonometric identities, we get This says that for a strongly attracting limit cycle, even if it is very close to the SNIC, the adjoint, and thus the PRC, can have a nearly sinusoidal shape. The validity of the approximation can be seen by looking at the relative sizes of the two terms that comprise

*z*in equation 4.2. The term

_{u}*z*multiplying is

_{r}*O*(1/

*a*), so that for

*a*large, the approximation will hold. However, as

*a*decreases, the term involving

*z*will become important. In fact, for any fixed

_{r}*a*(big or small), as , both

*z*and become large, but

_{r}*z*is much larger (

_{r}*z*is

_{r}*O*(1/(1−

*b*

^{2})) while ). Thus, the term with

*z*will eventually dominate no matter what value of

_{r}*a*is used. Indeed, as , and except near . Thus, independent of

*a*, as , the adjoint will look like , as the theory (Ermentrout, 1996; Brown et al., 2004; Canavier, 2006) predicts. The point is that the magnitude of

*z*is proportional to 1/[

_{r}*a*(1−

*b*

^{2})], and the magnitude of is proportional to so that for

*a*large, the

*z*term will dominate only for . In the remainder, we numerically explore the adjoint.

_{r}*s*is the amplitude of the perturbation in the

*u*direction. Comparison of Figure 3a, which shows the PRC for

*s*=0.1, and Figure 3b shows that the adjoint is still a fairly good approximation to the PRC even for a large perturbation of

*s*=0.1.

Having established that the adjoint is a reasonable approximation of the PRC, we plot the adjoints for different values of and *a*, for *b* near the SNIC value of 1. We have chosen *b* fairly close to the value of the saddle node. Despite this, we see in Figure 5, where , that for the strongly attracting limit cycle (*a*=100), the adjoint or PRC is nearly sinusoidal, in agreement with our asymptotic approximation. The approximation is predicated on the idea that the contribution of *z _{r}* is small for large

*a*. As our calculations above show, as

*a*decreases, this contribution increases and, depending on the value of , is seen through the appearance of a nearly nonnegative PRC or adjoint as predicted by the theory for systems near a saddle node.

With *a*=1, the negative region is a small part near 0 phase. The shape of the PRC curves in this case can be appreciated by the geometry of the isochrons in Figure 2. When *a*=1, the curvature of the isochrons in the neighborhood of leads to a shift to an isochron associated with a greater phase, and hence a positive PRC curve. However, as *a* increases, the isochrons have less curvature, and the PRC, as approximated by the adjoint, becomes more bimodal.

As with the case for *a*=1, for and strong attraction, *a*=10, the negative lobe of the adjoint will disappear as *b* gets closer and closer to *b*=1. Figure 5b illustrates this point. For *a*=100, we need before the asymptotic theory becomes valid.

## 5. Discussion

One of the compelling questions in the study of rhythms in the nervous system is how neurons synchronize. Weak coupling theory and other approximations have led to some general principles about what leads to neural synchronization (Van Vreeswijk, Abbott, & Ermentrout, 1994; Hansel et al., 1995; Pfeuty, Mato, Golomb, & Hansel, 2003; and Wang, 2010, for a recent review). One such principle is that nonnegative phase-resetting curves lead to synchrony with inhibitory synapses and are less likely to synchronize with excitatory synapses. Thus, one of the interesting biophysical questions is what determines the shape of the PRC. Excitatory synaptic coupling between neurons with a bimodal PRC encourages synchrony. Starting with the work in Hansel et al. (1995) and Ermentrout (1996), a common assumption has been that neurons that become rhythmic through a saddle node on a limit cycle (SNIC) will have a nonnegative PRC. This is based on theory that is strictly true only very close to the bifurcation when the stimulus will lead to a delay of the next action potential.

In this letter, we have described a very simple model that undergoes a SNIC bifurcation, and yet, quite close to the bifurcation, the PRC can be almost a perfect sine wave and thus strongly bimodal. If we regard the distance from the bifurcation as proportional to the drive, and thus the frequency of the oscillation, we see that the PRC of this model is very sensitive to a change in frequency (Gutkin, Ermentrout, & Reyes, 2005) and can switch from class I to class II over a very narrow range of parameters. Three parameters play a large role in determining the shape of the PRC; the attraction to the limit cycle *a*, the distance from the SNIC *b*, and finally the position of the stable equilibria for *b*>1, determined by . In this work, we find that for , the PRC will be largely positive, whereas for , the PRC will be largely negative. Consequently, the statement that systems with a SNIC will generally have a positive PRC is not true, since systems with a SNIC can also be bimodal or largely negative, as shown in Figures 3 and 5. However, it is nevertheless possible that real biological oscillators that have a largely positive PRC arise in general from an oscillation by a SNIC bifurcation.

Equations 3.2 are not suitable for detailed modeling of neural systems. However, the trace in Figure 4 is similar to the shape of an action potential. Further, the related model, in which there is uniform velocity on the unit circle (Hoppensteadt & Keener, 1982; Guevara & Glass, 1982), predicts many of the qualitative features observed when cardiac tissue is subjected to periodic stimulation. Similarly, we are optimistic that equations 3.2 might be useful in understanding the dynamics of nerve cells under periodic stimulation or when coupled to other neurons.

Our main conclusions can be summarized as follows. If the location of the incipient SNIC bifurcation is at a position of the limit cycle at which a stimulus will lead to a delay of the next action potential, the PRC will be primarily negative; if the location of the incipient SNIC bifurcation is at a position of the limit cycle at which a stimulus will lead to an advance of the next action potential, the PRC will be primarily positive; and if the location of the incipient SNIC bifurcation is at a position of the limit cycle lying between the region of advance and delay of the next action potential (e.g., for in Figure 3 here), the PRC will be primarily bimodal. However, in the last case, due to the curvature of the isochrons, the PRC may nevertheless be primarily positive in the asymptotic limit (see Figure 5).

## Acknowledgments

L. G. thanks the Natural Sciences Engineering and Research Council (Canada) and B. E. thanks the National Science Foundation (USA) for financial support.

## References

## Note

^{1}

Sometimes this is called type I rather than class I. Since type 1 and type 0 have also been used to refer to the topology of PRCs (Krogh-Madsen et al., 2012), we use the class to classify the shapes of the PRCs and type to classify the topology.