Table of Contents Author Guidelines Submit a Manuscript
ISRN Computational Biology
Volume 2013 (2013), Article ID 230571, 11 pages
Research Article

All Phase Resetting Curves Are Bimodal, but Some Are More Bimodal Than Others

Department of Physics and Astronomy, College of Charleston, 66 George Street, Charleston, SC 29424, USA

Received 22 August 2013; Accepted 19 November 2013

Academic Editors: P. Durrens and A. Fedorov

Copyright © 2013 Sorinel A. Oprisan. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Phase resetting curves (PRCs) are phenomenological and quantitative tools that tabulate the transient changes in the firing period of endogenous neural oscillators as a result of external stimuli, for example, presynaptic inputs. A brief current perturbation can produce either a delay (positive phase resetting) or an advance (negative phase resetting) of the subsequent spike, depending on the timing of the stimulus. We showed that any planar neural oscillator has two remarkable points, which we called neutral points, where brief current perturbations produce no phase resetting and where the PRC flips its sign. Since there are only two neutral points, all PRCs of planar neural oscillators are bimodal. The degree of bimodality of a PRC, that is, the ratio between the amplitudes of the delay and advance lobes of a PRC, can be smoothly adjusted when the bifurcation scenario leading to stable oscillatory behavior combines a saddle node of invariant circle (SNIC) and an Andronov-Hopf bifurcation (HB).

1. Introduction

Neural oscillators are excitable cells which means that as soon as a parameter, such as an external bias current or an ionic conductance, crosses a threshold, then the system switches from a stable rest state characterized by small dampened excursions of the membrane potential to a high amplitude excursion of the membrane potential, called action potential (AP). Repetitive firing excitable cells were categorized by Hodgkin [1] as Class 1, that is, the firing frequency is continuously tunable by adjusting a bias current and such neurons can fire with arbitrarily low frequencies, or Class 2, that is, discontinuous curve with a nonzero firing frequency threshold.

Excitable cells generate APs when a certain control parameter crosses a threshold; that is, a bifurcation occurs in their state space [2]. It became clear that threshold phenomena, such as the AP, are related to fundamental properties of the mathematical model describing the neuron [3] and its corresponding bifurcation diagram [4]. Class 1 excitability was ascribed to a saddle node of invariant circle (SNIC) bifurcation and Class 2 excitability to an Andronov-Hopf bifurcation (HB) (see [5, 6] and references therein). Since neurons are excitable cells operating near a bifurcation point, it is mathematically possible to significantly reduce the number of variables describing the state of a neuron to a few essential ones by using the so-called canonical forms (see [2] for an extensive review of existing canonical models in neuroscience). For example, the four-dimensional conductance based Hodgkin-Huxley (HH) model [1], which was originally introduced as a model of the giant squid axon, includes mathematical descriptions of voltage-gated inward (depolarizing) sodium and an outward (hyperpolarizing) potassium currents. However, using additional assumptions, for example, “instantaneous” activation/inactivation of some ionic channels, HH model can be reduced to two-dimensional models, such as FitzHugh [7] or Morris and Lecar [8]. In particular, Morris-Lecar model used in this paper relies on an instantaneously activating and voltage-gated inward current and a voltage-gated outward delayed rectifier channel (see the appendix for detailed equations). While dimension reduction results in a low-dimensional and sometimes even analytically tractable model, the qualitative aspects of neural dynamics are preserved in canonical models since near the bifurcation points some dimensions become trivial (see [6, 9] for a review of methods and applications of dimension reduction). The most significant benefit of dimension reduction is a better understanding of the general principles that govern AP generation through phase space analysis (see [4, 9] and references therein).

The characteristics of repetitive firing of neural oscillators, such as the intrinsic period , the firing threshold, or the duration of the AP, are determined by their biophysical details, such as the morphology and the ionic channels expressed. However, there are circumstances, for example, in the study of phase locked modes and synchrony of neural populations, when a phenomenological understanding of neural activity suffices. A phase resetting curve (PRC) is not concerned with the biophysical details of neurons but rather tabulates the transient changes in the firing period due to an external perturbation (usually a brief current or synaptic conductance pulse). The most significant effect of a perturbation occurs during the cycle that contains the perturbation and is quantified by the first order PRC; that is, , where is the intrinsic period of oscillation of isolated neuron and is the transiently modified duration of the current cycle due to a perturbation applied at the stimulus time or phase (see Figure 1(a)). Higher order PRCs are defined in a similar way in order to account for long-lasting effects of presynaptic inputs on the second and subsequent cycles (see [10] and the references therein).

Figure 1: The stable membrane potential oscillatory solution of Morris-Lecar model neuron (continuous line in panel (a)) has an intrinsic period , which can be transiently changed to (dashed line in panel (a)) due to a brief (inhibitory) perturbation applied at the stimulus time (vertical arrow) or stimulus phase . Equally, spaces phases, apart, are marked by solid circles along the unperturbed membrane potential trace (continuous line). The first order phase resetting curve (b), , represents the normalized change in the transient period of oscillation and is always bimodal. The amplitude, , of the delay lobe is very large compared to the amplitude, , of the advance lobe when the ML oscillator is very close to its SNIC bifurcation point (see inset in panel (b)).

Despite its limitations, the PRC is an invaluable tool in predicting the phase locked and synchronous modes of neural networks. The PRC-based approaches to treat Parkinson’s disease and epilepsy are at the forefront of practical applications of the extensive theory of phase resetting. It is believe that resting tremor is caused by a cluster of neurons located in the thalamus and the basal ganglia which fire synchronously [1113]. By using the PRC method, it is possible to predict the shape and timing of the necessary stimuli that destabilize a synchronous state of a large population of neurons with a single pulse [14], double pulses [15], or bipolar pulses [16]. Deep brain stimulation uses permanently implanted electrodes into a target area of the brain with the purpose of delivering precise current pulses aiming at correcting abnormal synchronization among neurons through careful phase resetting [1518]. Epileptic seizures can also be suppressed by applying brief current stimuli in order to reset the phases of some neurons and break the synchrony of the entire network [1921].

The shape of the PRC plays a crucial role in the way neurons synchronize or desynchronize when they are part of a neural network. It has been shown that neurons characterized by a Type II PRC, which has an almost sinusoidal shape and comparable sizes of the positive and negative lobes, synchronize more readily due to their ability to both slow down (if the oscillator is ahead) and speed up (if the oscillator is behind) [2224]. In contrast, the neurons endowed with a Type I PRC, which have a disproportionate ratio of the two PRC lobes, can predominantly speed up (or only slow down) their rhythm when coupled with other neurons. Figure 1(b) shows an example of a Type I PRC with a disproportionate ratio between the amplitudes of the positive and negative lobes. Traditionally, if the ratio of the amplitudes of the two lobes was less than [25] or other arbitrary number, then the PRC was classified as “unimodal” (see Figure 1(b)). The distinction between Type I and Type II PRCs is more profound than the simple visual appearance and was linked to the bifurcation mechanisms leading to stable oscillations.

It was long advocated that the PRCs are either of Type I (see Figure 1(b)) when the bifurcation mechanism leading to oscillatory behavior is a saddle node of invariant circle (SNIC) or of Type II when the bifurcation is of Andronov-Hopf (HB) type [2224, 26, 27]. Recent theoretical studies suggested that a SNIC bifurcation can also lead to bimodal PRCs [2830].

In this paper we shown that all PRCs are bimodal and the degree of bimodality of a PRC can be smoothly adjusted if the bifurcation scenario leading to stable limit cycle oscillations combines a SNIC and a HB. Furthermore, our current study extends previous results [28, 30] by introducing the concept of neutral points, that is, phases at which external perturbations produce no phase resetting, and showing that a SNIC that ends with a HB leads to a PRC that is a linear combination of both Type I and Type II PRCs.

2. Theory

2.1. Natural Coordinates for Limit Cycle Oscillators

Biological oscillators with a relatively stable periodic activity are mathematically described as limit cycle (LC) oscillators (see [28] and references therein). In order to gain insight into the oscillatory behavior of nonlinear oscillators and their corresponding PRCs, we used a natural reference frame, that is, the tangent and normal to the LC frame, instead of the usual fixed reference frame [29]. In such a natural reference frame, any external perturbation of the LC oscillator determines either a tangent displacement, that is, a sudden phase jump, or a normal (to the LC) perturbation. While the contribution of the tangent perturbations to the PRC is straightforward to account for, the phase resetting induced by normal perturbations is notoriously harder to account for [30, 31] and, according to [32], “the (A/N: canonical) model above may be of limited value since it fails to reflect the relaxation nature of the dynamics.”

Without reducing the generality of our results, here we only focused on planar neural oscillators whose vector flows are , where are the fixed reference frame variables, which in the case of Morris-Lecar model [27, 33] are the membrane potential, , and the potassium conductance, , and the superscript means transposed (see the appendix for model’s equations).

The analysis of the response of a dynamical system to external perturbations is significantly clearer in the moving reference frame of the figurative point traveling along the LC, , than in the fixed frame of the original system of differential equations (see Figure 2(a)). The choice of the unitary matrix that switches between the reference frames, that is, , is not unique [3436]. However, the unit vector of the tangent direction is uniquely defined by and then we chose one of the two possible orientations of the normal vector to be .

Figure 2: A stable limit cycle (LC) generates a closed trajectory in the fixed reference frame . (a) A membrane potential perturbation moves the figurative point from the unperturbed LC to . The tangent projection of the phase space perturbation is and represents an instantaneous geometric phase jump that corresponds to a temporal phase resetting. The normal projection of is and it relaxes back to the unperturbed limit cycle producing an additional geometric phase displacement that translates into a phase resetting. After a time , the figurative point would have traveled to along the unperturbed LC whereas the perturbed figurative point reached . The normal distance between the two end points is and the corresponding tangent displacement is the geometric phase resetting, [30], which corresponds to a (temporal) phase resetting measured by . (b) In the absence of any coupling between the normal and tangent directions of the LC, that is, , the only phase resetting is determined by . When the perturbation is perpendicular to the LC, it results in that and there is no phase resetting (points “” and “”) at these two neutral points.

In the fixed reference frame, , the temporal evolution of small perturbations, , from is described by the linearized equations of the dynamical system; that is, , where is the Jacobian matrix evaluated along . Using the unitary transform U, the same perturbation with respect to the moving reference frame, , is given by [29, 34] where , + , , and (we also corrected a typo in of [34]). The coefficient measures the effect along the direction for a perturbation applied along the direction with respect to . For example, the coefficient quantifies the effect of a pure normal perturbation measured along the tangent direction; that is, the amount of phase resetting induced along tangent direction to by relaxation processes. In all planar systems , which means that a perturbation tangent to never induces a normal displacement of the figurative point.

While the Jacobian matrix describes the stability of a solution to small perturbations with respect to the fixed reference frame , the matrix (see (1)) determines the stability with respect to normal and tangent perturbations at different phases along . A negative value of any coefficient means that the perturbation decays exponentially along the respective direction, that is, the LC is stable to perturbations applied at that phase.

In the PRC theory, the actual perturbation only affects the membrane potential equation; that is, with respect to the fixed reference frame, (see thick horizontal arrow in Figure 2(a)). The corresponding tangent and normal projections of the perturbation can be found using the unitary matrix : Equation (2) gives the initial conditions with which we integrated the perturbation equations (1). Since for any planar dynamical systems, , the differential equation for normal perturbations, that is, , has the analytic solution: where the initial condition is determined by the injected current of amplitude ; that is, . The differential equation for tangent perturbation, which actually determines the total amount of phase resetting due to the current perturbation , contains both pure tangent perturbations (through coefficient) and normal perturbations (through ); that is, which has the analytic solution: where .

Those fundamental theoretical results highlighted by (3) and (5) relate to the local stability coefficients of the normal and tangent geometrical displacements from the unperturbed limit cycle . In the subsequent section, we will confer these geometric (phase space) displacements into temporal advances or delays of the perturbed figurative point, that is, phase resetting.

3. Results

3.1. Why Are All PRCs Bimodal?

The fact that all PRCs are bimodal results from the existence of two remarkable points along —the neutral points or zero phase resetting points. A perturbation applied at a neutral point induces no phase resetting. For planar systems, the existence of the two neutral points can be explain in terms of the stability coefficients of the LC. Indeed, let us assume for a moment that the crossover coefficient vanishes. This means that the only contribution to the PRC is due to tangent perturbations through . When the applied perturbation is perpendicular to , its corresponding tangent projection vanishes; therefore, its contribution to the PRC also vanishes because the initial tangent displacement in (5) is zero. Since the normal perturbation does not contribute to the tangent displacement (because ), the points where the brief (infinitesimal) voltage perturbation is normal to are the neutral points, that is, the points “” and “” in Figure 2(b). In other words, if , then the maximum and the minimum of the membrane potential are the neutral points of the limit cycle.

If , the exact positions (or phases) of the two neutral points and in Figure 2(b) can be determined from the requirement that the contribution of the pure tangential perturbations to (5), that is, , is exactly balanced by the contribution of the normal relaxation processes through . Furthermore, for planar systems, there are two and only two neutral points along ; that is, the PRC is only bimodal. Indeed, assuming again for a moment that there is no normal-to-tangent coupling between perturbations (), then the only contribution to the PRC is due to . On the other hand, contribution vanishes at local minima or maxima of the membrane potential because the tangent projection of the perturbation, , vanishes. For a planar system, there is only one maximum and only one minimum of the membrane potential; otherwise at least one more independent variable (in addition to the membrane potential and one slow variable ) is required in order to further modulate the membrane potential and produce more local minima and maxima, such as bursting activity. Furthermore, at neutral points the sign of the PRC flips. Indeed, assuming for a moment no normal-to-tangent perturbation coupling (), then the tangent projection of a membrane potential perturbation (see Figure 2(b)) always delays the next spike if applied during the repolarization of the AP around the two neutral points of the limit cycle (see point in Figure 2(b)). On the other hand, outside the repolarization region of the AP, a brief depolarization always advances the next spike. Therefore, the PRC must flip its sign at neutral points.

3.2. Are PRCs Unimodal Near a SNIC Bifurcation?

The theoretical arguments of our conjecture regarding the bimodality of all PRCs are supported by numerical simulations. We found that a brief depolarization applied during the repolarization phase of the AP always delays the next spike of ML model neuron (lengthening the duration of the current cycle), whereas a brief hyperpolarization during the repolarization phase always advances the next spike (see Figure 3). All PRCs in Figure 3 were numerically obtained by integrating the differential equations of the ML model (see the appendix) with a current perturbation in the shape of a rectangular pulse of amplitude (in absolute units, which is 2.45% of the bias current at the SNIC bifurcation point) and a duration of .

Figure 3: All PRCs, including Type I, are bimodal. The bimodality increases while the system moves away from the SNIC bifurcation point that generated the stable limit cycle oscillations. Very close to the SNIC bifurcation point (A1) , the positive portion (delay of the next spike induced by a depolarization) is not even visible, unless we magnify it 1000 times (solid dots in panel (B1)). The bimodality of the PRC becomes more transparent away from the SNIC bifurcation; for example, (B1), (C1), and (D1), respectively. The amplitude of the positive lobe of the PRC is magnified 1000 times (large solid circles) and plotted together with the action potential (dashed line). The two neutral points where no resetting occurs are both on the repolarizing slope of the action potential (see the vertical dashed lines in all panels of the second column).

Close to the SNIC bifurcation point (see Figure 3(A1)) the positive portion of the PRC is hardly visible in print unless it is magnified 1000 times (see the solid dots in Figure 3(A2)). When the magnified PRC is plotted together with the AP (Figure 3(A2)), we notice that the delay is only induced during a small portion of the repolarization phase that starts shortly after the membrane potential maximum (the first neutral point in the absence of normal-to-tangent coupling) and ends before the membrane potential minimum (the second neutral point in the absence of normal-to-tangent coupling). We also notice from Figure 3 that, very close to the SNIC bifurcation point, small perturbations induce significant phase resettings. Indeed, stimulus amplitude has a very large effect (over maximum phase resetting) close to the SNIC bifurcation point (Figure 3(A1)), whereas for a stimulus amplitude that is only about away from the SNIC bifurcation the amount of phase resetting is reduced 100 times (Figure 3(A2)). The first column of Figure 3 clearly supports the conjecture that Type I PRC, which is characteristic to stable oscillations emerging through a SNIC bifurcation, is always bimodal and its bimodality becomes more apparent while the system moves away from the SNIC bifurcation point. At the same time, the positive portion of the PRC, that is, the delay of the subsequent spike in response to a brief depolarization, always occurs during the repolarization phase (second column in Figure 3), as we conjectured.

The curve shows the characteristic Class I excitability behavior (Figure 4(a)); that is, the firing frequency can be continuously adjusted arbitrarily close to zero [1], when the bias current approaches the SNIC bifurcation, up to about 120 Hz (Figure 4(a)). We used the freely available software package AUTO [37] to obtain the bifurcation diagram (Figure 4(b)) for a fixed value of calcium conductance () and variable bias current. The bifurcation diagram has a branch of hyperpolarized stable steady states for low values of bias current (continuous line), which becomes unstable at the SNIC bifurcation point (see point in Figure 4(b)). The emerging LC at SNIC bifurcation is stable with the peak-to-peak amplitude represented by the pairs of solid circles for each value of the bias current up to the maximum that corresponds to the limit point in Figure 4(b). To further support our conjecture that all PRCs are bimodal, we numerically computed the ratio of the PRC lobes (see Figure 1(b)). Very close to the SNIC bifurcation point the amplitude of the smaller lobe is a few orders of magnitude smaller than the magnitude of the other lobe (Figure 4(c)), which justifies the name “unimodal” for Type I PRC, associated with SNIC bifurcations.

Figure 4: The firing frequency of Type I excitable cells can be continuously adjusted over a wide range by changing the bias current and can be tuned arbitrarily close to zero (a). The bifurcation diagram for a fixed calcium conductance (b) shows a clear SNIC bifurcation with the fold characterized by the two turning points and , which occur at the SNIC bifurcation point. In all bifurcation diagrams, solid lines mark stable steady states and dashed lines mark unstable ones. Symmetric pairs of solid circles for the same bias current mark peak-to-peak voltages of the stable oscillations, while the open circles mark unstable limit cycle oscillations. On the depolarized steady state branch (continuous line for in panel (b)) there is a subcritical Andronov-Hopf and double-limit cycle bifurcations with stable oscillations emerging at (see the point ) while the corresponding steady state is still stable. The loss of stability is due to a double-limit cycle bifurcation, characterized by the simultaneous appearance of two limit cycles of opposite stability (solid circles—stable and open circles—unstable limit cycles). The local Andronov-Hopf bifurcation, also named degenerated Andronov-Hopf bifurcation [38, 39], occurs at (HP point in panel (b)). The ratio of the amplitudes of the two lobes of the PRC can be continuously tuned over a few orders of magnitude (c) from almost (near the SNIC bifurcation point) up to almost one unit (close to Andronov-Hopf bifurcation marked as HB). The numerically generated PRCs in response to brief square pulses were fitted with a linear combination of the theoretical Type I PRC () and Type II PRC (). The ratio of the coefficient of against shows that very close to the SNIC bifurcation point the Type I PRC “unimodal” character dominates over the Type II “bimodal” character of the PRC (d). When the ratio of the two parameters is close to 1 (see the horizontal dashed line in panel (d)), the two types balance. The combination of a SNIC bifurcation point at low values of the bias current and Andronov-Hopf bifurcation point at large bias currents (SNIC-HB bifurcation diagram) leads to a PRC that is a linear combination of a pure Type I and a pure Type II PRC.
3.3. Why the Bimodality of the PRCs Can Be Smoothly Tuned?

While it is clear that every PRC is bimodal, one question is still open: why does the bimodality of the PRCs change with the control parameters? For example, in Figure 4(c) we notice that the ratio of the amplitudes of the two lobes smoothly increases with the bias current until they are comparable.

How is it possible for the PRC of a system for which the oscillatory behavior emerges through a SNIC bifurcation (and whose PRC was supposed to be “unimodal”) to be so bimodal that we cannot distinguish it from a Type II PRC (which is characteristic for HB)? We believe that the answer must take into account the entire structure of the bifurcation diagram and not only the fact that at some point the stable oscillatory behavior was generated by a SNIC. For example, the complete bifurcation diagram for ML model shown in Figure 4(b) has two bifurcation points, both of which generate stable limit cycle oscillations: on the hyperpolarized stable steady state branch there is indeed a SNIC bifurcation point for low values of the bias current, and on the depolarized branch there is an Andronov-Hopf (HB) point (Figure 4(b)). The oscillatory behavior emerges through a SNIC bifurcation by slowly increasing the bias current along the hyperpolarized branch ( point in Figure 4(b)), which produces a PRC that can be approximated with the “unimodal” Type I expression in the close proximity of the SNIC bifurcation [2, 22, 26]. On the other hand, stable oscillations emerge as soon as the HB point is reached (HB point in Figure 4(b)) by slowly decreasing the bias current along the depolarized branch. As such, the close proximity of the HB produces a PRC that can be approximated with the “bimodal” Type II expression with being a constant phase [2, 22, 26]. We conjectured that no PRC is purely of Type I or Type II but a (linear) combination of the two types. To check this conjecture, we fitted the numerically generated PRCs with a linear combination of the and ; that is, , where and are fitting coefficients. Ideally, we expect that the PRC is entirely captured by a Type I iPRC very close to a SNIC bifurcation; that is, .

By the same token, we expect that the PRC is entirely captured by a Type II iPRC very close to an Andronov-Hopf bifurcation point; that is, . Since the amplitudes and depend on the characteristics of the injected current perturbation (amplitude and duration), we computed their ratio , which is insensitive to such variations (see Figure 4(d)). Our numerical simulations show that indeed the contribution of is one order of magnitude larger than that of near a SNIC bifurcation point whereas the situation is reversed near a HB point (see Figure 4(d)). Furthermore, by setting an arbitrary ratio for classifying the PRCs into “unimodal” or “bimodal”, for example, 0.175 [25] or , we can actually use Figure 4(d) to estimate the range of the bias currents for which such an arbitrary classification is valid.

3.4. How Robust Is the SNIC-HB Scenario That Produces Tunable PRC Bimodality?

The complete bifurcation diagram shown in Figure 4(b) was obtained for fixed calcium conductance (). Since the pair SNIC-HB of bifurcation points is crucial to a tunable bimodality of the PRC, we also explored the reproducibility of such a bifurcation diagram while changing . For example, for large values of calcium conductance, , the bifurcation scenario is always the same: while increasing the bias current from very low values (below 0.08), stable LC oscillations emerge through a SNIC bifurcation (point in Figure 5(a)), but while decreasing the bias current from large values (above 0.3), stable limit cycle oscillations emerge through an Andronov-Hopf bifurcation (HB point in Figure 5(a)). It is this particular bifurcation scenario, with a SNIC bifurcation at one end and an Andronov-Hopf bifurcation at the other end of the bifurcation diagram, that leads to a PRC which is a linear combination of a Type I iPRC (typically near a SNIC bifurcation) and a Type II iPRC (typically near an Andronov-Hopf bifurcation) as shown above. For a critical value of calcium conductance (close to 0.6), both ends of the stable LC oscillations are marked by Andronov-Hopf bifurcation points (Figure 5(b)). The SNIC-HB bifurcation scenario changes to a HB-HB when the two limit points and of the SNIC fold collide and give way to a HB point (see Figure 5(b)). The bifurcation diagrams (Figure 5) also show that large calcium conductances lead to larger peak-to-peak membrane potential oscillations, especially through a significant increase of the depolarization overshoot. This effect is due to the fact the a larger calcium conductance drives the depolarization of the cell more effectively as soon as the membrane potential exceeds the half activation threshold of calcium channels (the parameter of the model described in the appendix) and leads to a rapid and strong overshoot (Figure 5(a)). Large calcium currents, due to a large value of , require a relatively modest depolarizing bias current compensation to sustain oscillations. This is the reason why the bias current range is narrower for larger in the SNIC-HB bifurcation diagrams. On the other hand, for small calcium conductances (Figure 5(b)) the oscillations have small amplitude since calcium channels do not reach their half activation voltage and the dynamics of the system are dominated by potassium channels. Indeed, small hyperpolarizations induced by potassium immediately and dramatically decrease the number of open voltage-gated calcium channels (since the calcium channels already work below their half activation threshold). At the same time, the bias current required to compensate the dynamics described by a HB-HB scenario, which is dominated by the outward potassium current, is significantly larger than the bias required for a SNIC-HB scenario (compare bias current values in Figure 5(a) with those in Figure 5(b)).

Figure 5: Bifurcation diagrams of ML model neuron for variable bias current at high (a) and low (b) calcium conductances, respectively. For large , stable limit cycle oscillations emerge at low bias currents through a SNIC bifurcation (point ) and through Andronov-Hopf (HB point) for large bias currents (a). The stable oscillations that emerge in a SNIC-HB diagram have a PRC that is a linear combination of a pure Type I PRC (characteristic to SNIC bifurcation) and pure Type II PRC (characteristic to HB). At smaller calcium conductance, the two limit points and of the SNIC bifurcation collapse into a HB point and a HB-HB diagram emerges (b). The peak-to-peak amplitude of oscillations in a HB-HB diagram is smaller than that in a SNIC-HB diagram due to a reduce activation of calcium channels and the fact that the dynamics of the system are dominated by the outward potassium current. As a result, in a HB-HB system the bias current is larger than that in the case of a SNIC-HB system to compensate for the larger outward potassium current.

4. Discussions

Network models showed that selectively blocking some channels in critical inhibitory neurons can lead to switches from unsynchronized high frequency to synchronized low frequency network oscillations [40]. General anesthetics are selective blockers that could also induce dramatic change in the PRC shape and switch the brain from firing patterns associated with consciousness to states associated with general anesthesia [4143]. Therefore, the conductances of voltage-gated channels, such as calcium in our ML example, could be pharmacologically manipulated to produce a dramatic change in the bimodality of the PRC. Although this paper focuses on understanding the structure of PRCs and the effect of biophysically relevant control parameters, such as the bias current and calcium conductance, at the single cell level, the results are relevant actually at network level. In the end, the PRC is a tool used for predicting neural population behavior [14, 44].

We analyzed the mechanisms leading to a smooth and continuous change of PRC’s bimodality in a planar model. We used a ML model neuron not only because it represents a biophysically realistic model of an actual neuron [8, 27, 33], but also because, by changing only two parameters, that is, the calcium conductance and the bias current, it displays both SNIC and Andronov-Hopf bifurcations. Some authors used the ratio of the positive and the negative amplitudes of the PRCs to classify them as unimodal or bimodal (see [26, 28] and references therein). In this paper, we showed that all PRCs are bimodal regardless of the bifurcation mechanism that led to oscillatory behavior. Moreover, we showed that there are two remarkable points on any planar limit cycle, which we called the neutral or zero resetting points, at which the PRC flips its sign (see Figure 2). Although the phases of the two neutral points slightly change with the parameters of the external perturbation (amplitude and duration of stimulus), their location can be inferred from the model’s equations using the stability coefficients of the LC. The phase at which a brief depolarization can induce a delay (positive resetting) is situated along the repolarization region of the AP (see Figure 3). Since for fast spiking neurons the repolarization covers only a minute fraction of the [30], it may raise concerns regarding real-world applicability of such weak bimodality. Experimental studies showed that a brief pulse can desynchronize a fully synchronized cluster of neural oscillators if the pulse hits the cluster in a very narrow phase range (5% ) [14, 16, 45]. It is also possible for the two neutral points to actually collide, in which case the PRC becomes truly unimodal. We showed that a ML model can mimic Class I excitability (with a continuous frequency versus current curve—see Figure 4(a)) with a SNIC-HB bifurcation diagram (Figure 4(a)). The tunable bimodality of the PRC is a result of the presence of both SNIC and HB, which allows fine mixing of Type I character (specific to a SNIC bifurcation) with Type II character (specific to a HB) into the final PRC (Figure 4(c)).

Oscillatory activity is a hallmark of neuronal network function in various brain regions, including the olfactory bulb [4648], thalamus [49, 50], hippocampus [51, 52], and neocortex [53]. The existence and stability of synchronous and phase locked modes can be predicted using PRCs (see [10, 31, 5456] and references therein). Neural networks oscillate over more than three orders of magnitude and include delta (0.5–3 Hz), theta (3–8 Hz), gamma (30–90 Hz), and ultrafast (90–200 Hz) oscillations [57, 58]. Stable oscillatory activity at single cell or population levels was associated with higher brain functions [59, 60], such as gamma rhythm that is believed to be the basis of temporal encoding [61, 62], sensory binding of features into a coherent perception [59, 63], and storage and recall of information [6466]. Disruption of synchronous rhythms underlies some psychiatric disorders, such as schizophrenia [67, 68]. The ability to smoothly change the PRC from a Type I to a Type II through experimentally available control parameters, such as bias currents or pharmacologically controlled ionic channels, has potential application in shaping the activity of neural cells and networks.


We used a dimensionless Morris-Lecar model [27, 33] described by the following equations: where is the membrane potential, is the slow potassium activation, and all ionic currents are described by , where is the conductance of the voltage-gated channel and is the corresponding reversal potential. In particular, the calcium current is , the potassium current is , and the leak current is . The reversal potentials for calcium, potassium, and leak currents are , , and , respectively. The steady state activation function for calcium channels is , where , , the steady state activation function for potassium channels is , where , , the inverse time constant of potassium channels is , the potassium and leak conductances are , , respectively, and the .

The two control parameters that can switch the bifurcation diagram from a SNIC-HB type to a HB-HB type (see Figure 5) are the calcium conductance and the bias current . For example, if and , (A.1) describe what was classified by Hodgkin as Class I excitable cells. For example, if and , (A.1) describe what was classified by Hodgkin as Class II excitable cells.

The two nullclines of (A.1) are two-dimensional curves and , respectively. A fixed point of (A.1) is any crossing of the two nullclines, that is, phase space points where both the rate of change of membrane potential and the potassium activation vanish. The linear stability of the fixed points is determined by the two eigenvalues of the Jacobian matrix of (A.1). Stable steady states have eigenvalues with negative real part and are represented by a solid line in any bifurcation diagram (see Figures 4(b) and 5). Saddle or hyperbolic points have real eigenvalues with opposed signs, whereas Andronov-Hopf points have pairs of purely imaginary eigenvalues.


This research was supported by the National Science Foundation Grant CAREER Award IOS 1054914.


  1. A. Hodgkin, “The local electric changes associated with repetitive action in a non-medullated axon,” The Journal of Physiology, vol. 107, pp. 165–181, 1948. View at Google Scholar
  2. E. M. Izhikevich, Systems in Neuroscience: The Geometry of Excitability and Bursting, The MIT Press, 2007.
  3. R. FitzHugh, “Mathematical models of threshold phenomena in the nerve membrane,” The Bulletin of Mathematical Biophysics, vol. 17, no. 4, pp. 257–278, 1955. View at Publisher · View at Google Scholar · View at Scopus
  4. J. Rinzel and G. Ermentrout, Methods in Neuronal Modeling, Chapter Analysis of Neural Excitability and Oscillations, The MIT Press, Cambridge, UK, 1989.
  5. S. A. Prescott, Y. De Koninck, and T. J. Sejnowski, “Biophysical basis for three distinct dynamical mechanisms of action potential initiation,” PLoS Computational Biology, vol. 4, no. 10, Article ID e1000198, 2008. View at Publisher · View at Google Scholar · View at Scopus
  6. R. M. Smeal, G. B. Ermentrout, and J. A. White, “Phase-response curves and synchronized neural networks,” Philosophical Transactions of the Royal Society B, vol. 365, no. 1551, pp. 2407–2422, 2010. View at Publisher · View at Google Scholar · View at Scopus
  7. R. FitzHugh, “Impulses and physiological states in theoretical models of nerve membrane,” Biophysical Journal, vol. 1, pp. 445–466, 1961. View at Google Scholar
  8. C. Morris and H. Lecar, “Voltage oscillations in the barnacle giant muscle fiber,” Biophysical Journal, vol. 35, no. 1, pp. 193–213, 1981. View at Google Scholar · View at Scopus
  9. F. C. Hoppensteadt and E. M. Izhikevich, Weakly Connected Neural Networks, Springer, New York, NY, USA, 1997.
  10. S. Oprisan, A Geometric Approach to Phase Resetting Estimation Based on Mapping Temporal to Geometric Phase, vol. 6 of Springer Series in Computational Neuroscience, 2012.
  11. R. Llinás and H. Jahnsen, “Electrophysiology of mammalian thalamic neurones in vitro,” Nature, vol. 297, no. 5865, pp. 406–408, 1982. View at Publisher · View at Google Scholar · View at Scopus
  12. J. Sarnthein and D. Jeanmonod, “High thalamocortical theta coherence in patients with Parkinson's disease,” Journal of Neuroscience, vol. 27, no. 1, pp. 124–131, 2007. View at Publisher · View at Google Scholar · View at Scopus
  13. A. Schnitzler and J. Gross, “Normal and pathological oscillatory communication in the brain,” Nature Reviews Neuroscience, vol. 6, no. 4, pp. 285–296, 2005. View at Publisher · View at Google Scholar · View at Scopus
  14. P. A. Tass, Phase Resetting in Medicine and Biology-Stochastic Modelling and Data Analysis, Springer, 1999.
  15. P. A. Tass, “Desynchronizing double-pulse phase resetting and application to deep brain stimulation,” Biological Cybernetics, vol. 85, no. 5, pp. 343–354, 2001. View at Publisher · View at Google Scholar · View at Scopus
  16. P. A. Tass, “Desynchronization of brain rhythms with soft phase-resetting techniques,” Biological Cybernetics, vol. 87, no. 2, pp. 102–115, 2002. View at Publisher · View at Google Scholar · View at Scopus
  17. A. L. Benabid, P. Pollak, C. Gervason et al., “Long-term suppression of tremor by chronic stimulation of the ventral intermediate thalamic nucleus,” The Lancet, vol. 337, no. 8738, pp. 403–406, 1991. View at Publisher · View at Google Scholar · View at Scopus
  18. P. A. Tass, “A model of desynchronizing deep brain stimulation with a demand-controlled coordinated reset of neural subpopulations,” Biological Cybernetics, vol. 89, no. 2, pp. 81–88, 2003. View at Publisher · View at Google Scholar · View at Scopus
  19. J. Engel, T. Pedley, and J. Aicardi, “Epilepsy: a comprehensive textbook,” in Epilepsy: A Comprehensive Textbook, vol. 1, Wolters Kluwer Health/Lippincott Williams & Wilkins, 2008. View at Google Scholar
  20. O. V. Popovych, C. Hauptmann, and P. A. Tass, “Control of neuronal synchrony by nonlinear delayed feedback,” Biological Cybernetics, vol. 95, no. 1, pp. 69–85, 2006. View at Publisher · View at Google Scholar · View at Scopus
  21. M. G. Rosenblum and A. S. Pikovsky, “Controlling synchronization in an ensemble of globally coupled oscillators,” Physical Review Letters, vol. 92, no. 11, Article ID 114102, 4 pages, 2004. View at Publisher · View at Google Scholar · View at Scopus
  22. B. Ermentrout, M. Pascal, and B. Gutkin, “The effects of spike frequency adaptation and negative feedback on the synchronization of neural oscillators,” Neural Computation, vol. 13, no. 6, pp. 1285–1310, 2001. View at Publisher · View at Google Scholar · View at Scopus
  23. B. S. Gutkin, G. B. Ermentrout, and A. D. Reyes, “Phase-response curves give the responses of neurons to transient inputs,” Journal of Neurophysiology, vol. 94, no. 2, pp. 1623–1635, 2005. View at Publisher · View at Google Scholar · View at Scopus
  24. D. Hansel, G. Mato, and C. Meunier, “Synchrony in excitatory neural networks,” Neural computation, vol. 7, no. 2, pp. 307–337, 1995. View at Google Scholar · View at Scopus
  25. T. Tateno and H. P. C. Robinson, “Phase resetting curves and oscillatory stability in interneurons of rat somatosensory cortex,” Biophysical Journal, vol. 92, no. 2, pp. 683–695, 2007. View at Publisher · View at Google Scholar · View at Scopus
  26. E. Brown, J. Moehlis, and P. Holmes, “On the phase reduction and response dynamics of neural oscillator populations,” Neural Computation, vol. 16, no. 4, pp. 673–715, 2004. View at Publisher · View at Google Scholar · View at Scopus
  27. B. Ermentrout, “Type I membranes, phase resetting curves, and synchrony,” Neural Computation, vol. 8, no. 5, pp. 979–1001, 1996. View at Google Scholar · View at Scopus
  28. G. B. Ermentrout, L. Glass, and B. E. Oldeman, “The shape of phase-resetting curves in oscillators with a saddle node on an invariant circle bifurcation,” Neural Computation, vol. 24, pp. 3111–3125, 2012. View at Google Scholar
  29. S. A. Oprisan, “Local linear approximation of the jacobian matrix better captures phase resetting of neural limit cycle oscillators,” Neural Computation, vol. 26, no. 1, pp. 1–26, 2014. View at Google Scholar
  30. S. A. Oprisan and C. C. Canavier, “The influence of limit cycle topology on the phase resetting curve,” Neural Computation, vol. 14, no. 5, pp. 1027–1057, 2002. View at Publisher · View at Google Scholar · View at Scopus
  31. S. A. Oprisan, A. A. Prinz, and C. C. Canavier, “Phase resetting and phase locking in hybrid circuits of one model and one biological neuron,” Biophysical Journal, vol. 87, no. 4, pp. 2283–2298, 2004. View at Publisher · View at Google Scholar · View at Scopus
  32. E. M. Izhikevich, “Neural excitability, spiking and bursting,” International Journal of Bifurcation and Chaos in Applied Sciences and Engineering, vol. 10, no. 6, pp. 1171–1266, 2000. View at Google Scholar · View at Scopus
  33. H. Lecar and R. Nossal, “Theory of threshold fluctuations in nerves. II. Analysis of various sources of membrane noise,” Biophysical Journal, vol. 11, no. 12, pp. 1068–1084, 1971. View at Google Scholar · View at Scopus
  34. F. Ali and M. Menzinger, “On the local stability of limit cycles,” Chaos, vol. 9, no. 2, pp. S1054–S1500, 1999. View at Google Scholar
  35. R. Bellman, introduction to Matrix Analysis. Society for industrial and Applied Mathematics, University City Service Center, Philadelphia, Pa, USA, 1997.
  36. T. S. Shores, Applied Linear Algebra and Matrix Analysis, Springer, 2007.
  37. B. Ermentrout, Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT for Researchers and Students, vol. 14, Society for Industrial and Applied Mathematics, Philadelphia, Pa, USA, 2002.
  38. J. Guckenheimer and I. S. Labouriau, “Bifurcation of the Hodgkin and Huxley equations: a new twist,” Bulletin of Mathematical Biology, vol. 55, no. 5, pp. 937–952, 1993. View at Publisher · View at Google Scholar · View at Scopus
  39. F. Takens, “Forced oscillations and bifurcations,” Communications of the Mathematical Institute, Rijksuniversiteit Utrecht, vol. 3, pp. 1–59, 1974. View at Google Scholar
  40. E. M. Izhikevich, “Subcritical elliptic bursting of Bautin type,” SIAM Journal on Applied Mathematics, vol. 60, no. 2, pp. 503–535, 2000. View at Google Scholar · View at Scopus
  41. P. Århem and C. Blomberg, “Ion channel density and threshold dynamics of repetitive firing in a cortical neuron model,” BioSystems, vol. 89, no. 1–3, pp. 117–125, 2007. View at Publisher · View at Google Scholar · View at Scopus
  42. P. Århem, G. Klement, and C. Blomberg, “Channel density regulation of firing patterns in a cortical neuron model,” Biophysical Journal, vol. 90, no. 12, pp. 4392–4404, 2006. View at Publisher · View at Google Scholar · View at Scopus
  43. H. Zeberg, C. Blomberg, and P. Århem, “Ion channel density regulates switches between regular and fast spiking in soma but not in axons,” PLoS Computational Biology, vol. 6, no. 4, Article ID e1000753, 2010. View at Publisher · View at Google Scholar · View at Scopus
  44. A. T. Winfree, The Geometry of Biological Time, Springer, 2001.
  45. P. A. Tass, “Stochastic phase resetting: a theory for deep brain stimulation,” Progress of Theoretical Physics Supplement, no. 139, pp. 301–313, 2000. View at Google Scholar · View at Scopus
  46. J. Beshel, N. Kopell, and L. M. Kay, “Olfactory bulb gamma oscillations are enhanced with task demands,” Journal of Neuroscience, vol. 27, no. 31, pp. 8358–8365, 2007. View at Publisher · View at Google Scholar · View at Scopus
  47. L. M. Kay and P. Lazzara, “How global are olfactory bulb oscillations?” Journal of Neurophysiology, vol. 104, no. 3, pp. 1768–1773, 2010. View at Publisher · View at Google Scholar · View at Scopus
  48. C. Martin, D. Houitte, M. Guillermier, F. Petit, G. Bonvento, and H. Gurden, “Alteration of sensory-evoked metabolic and oscillatory activities in the olfactory bulb of glast-deficient mice,” Frontiers in Neural Circuits, vol. 6, article 1, 2012. View at Google Scholar · View at Scopus
  49. S. W. Hughes, M. L. Lorincz, H. R. Parri, and V. Crunelli, “Infraslow (< 0.1Hz) oscillations in thalamic relay nuclei. Basic mechanisms and significance to health and disease states,” Progress in Brain Research, vol. 193, pp. 145–162, 2011. View at Publisher · View at Google Scholar · View at Scopus
  50. M. Steriade, E. G. Jones, and R. R. Llinas, Thalamic Oscillations and Signaling, John Wiley and Sons, Oxford, UK, 1990.
  51. G. Buzsáki, “Theta oscillations in the hippocampus,” Neuron, vol. 33, no. 3, pp. 325–340, 2002. View at Publisher · View at Google Scholar · View at Scopus
  52. M. Pignatelli, A. Beyeler, and X. Leinekugel, “Neural circuits underlying the generation of theta oscillations,” Journal of Physiology Paris, vol. 106, no. 3-4, pp. 81–92, 2011. View at Publisher · View at Google Scholar · View at Scopus
  53. G. A. Worrell, L. Parish, S. D. Cranstoun, R. Jonas, G. Baltuch, and B. Litt, “High-frequency oscillations and seizure generation in neocortical epilepsy,” Brain, vol. 127, no. 7, pp. 1496–1506, 2004. View at Publisher · View at Google Scholar · View at Scopus
  54. S. A. Oprisan, “Stability of synchronous oscillations in a periodic network,” International Journal of Neuroscience, vol. 119, no. 4, pp. 482–491, 2009. View at Publisher · View at Google Scholar · View at Scopus
  55. S. A. Oprisan, “Existence and stability criteria for phase-locked modes in ring neural networks based on the spike time resetting curve method,” Journal of Theoretical Biology, vol. 262, no. 2, pp. 232–244, 2010. View at Publisher · View at Google Scholar · View at Scopus
  56. S. A. Oprisan, “Existence and stability criteria for phase-locked modes in ring neural networks based on the spike time resetting curve method,” Journal of Theoretical Biology, vol. 262, no. 2, pp. 232–244, 2010. View at Publisher · View at Google Scholar · View at Scopus
  57. M. Bartos, I. Vida, and P. Jonas, “Synaptic mechanisms of synchronized gamma oscillations in inhibitory interneuron networks,” Nature Reviews Neuroscience, vol. 8, no. 1, pp. 45–56, 2007. View at Publisher · View at Google Scholar · View at Scopus
  58. G. Buzsáki and A. Draguhn, “Neuronal olscillations in cortical networks,” Science, vol. 304, no. 5679, pp. 1926–1929, 2004. View at Publisher · View at Google Scholar · View at Scopus
  59. C. M. Gray and W. Singer, “Stimulus-specific neuronal oscillations in orientation columns of cat visual cortex,” Proceedings of the National Academy of Sciences of the United States of America, vol. 86, no. 5, pp. 1698–1702, 1989. View at Google Scholar · View at Scopus
  60. U. Ribary, A. A. Ioannides, K. D. Singh et al., “Magnetic field tomography of coherent thalamocortical 40-Hz oscillations in humans,” Proceedings of the National Academy of Sciences of the United States of America, vol. 88, no. 24, pp. 11037–11041, 1991. View at Google Scholar · View at Scopus
  61. G. Buzsaki and J. J. Chrobak, “Temporal structure in spatially organized neuronal ensembles: a role for interneuronal networks,” Current Opinion in Neurobiology, vol. 5, no. 4, pp. 504–510, 1995. View at Publisher · View at Google Scholar · View at Scopus
  62. J. J. Hopfield, “Pattern recognition computation using action potential timing for stimulus representation,” Nature, vol. 376, no. 6535, pp. 33–36, 1995. View at Google Scholar · View at Scopus
  63. C. E. Schroeder and P. Lakatos, “Low-frequency neuronal oscillations as instruments of sensory selection,” Trends in Neurosciences, vol. 32, no. 1, pp. 9–18, 2009. View at Publisher · View at Google Scholar · View at Scopus
  64. J. K. Kleen, E. X. Wu, G. L. Holmes, R. C. Scott, and P.-P. Lenck-Santini, “Enhanced oscillatory activity in the hippocampal-prefrontal network is related to short-term memory function after early-life seizures,” Journal of Neuroscience, vol. 31, no. 43, pp. 15397–15406, 2011. View at Publisher · View at Google Scholar · View at Scopus
  65. J. E. Lisman, “Relating hippocampal circuitry to function: recall of memory sequences by reciprocal dentate-CA3 interactions,” Neuron, vol. 22, no. 2, pp. 233–242, 1999. View at Google Scholar · View at Scopus
  66. J. E. Lisman and M. A. P. Idiart, “Storage of 7 ± 2 short-term memories in oscillatory subcycles,” Science, vol. 267, no. 5203, pp. 1512–1515, 1995. View at Google Scholar · View at Scopus
  67. D. A. Lewis, T. Hashimoto, and D. W. Volk, “Cortical inhibitory neurons and schizophrenia,” Nature Reviews Neuroscience, vol. 6, no. 4, pp. 312–324, 2005. View at Publisher · View at Google Scholar · View at Scopus
  68. K. M. Spencer, P. G. Nestor, M. A. Niznikiewicz, D. F. Salisbury, M. E. Shenton, and R. W. McCarley, “Abnormal neural synchrony in schizophrenia,” Journal of Neuroscience, vol. 23, no. 19, pp. 7407–7411, 2003. View at Google Scholar · View at Scopus