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.
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.
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 Tj−Tj−1, j>1 differs from T depends on the rate at which the rhythm returns to the limit cycle oscillation.
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.
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 .
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.
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, .
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 zr 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.
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).
L. G. thanks the Natural Sciences Engineering and Research Council (Canada) and B. E. thanks the National Science Foundation (USA) for financial support.
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.