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

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 timingofthestimulus.Weshowedthatanyplanarneuraloscillatorhastworemarkablepoints,whichwecalledneutralpoints,wherebriefcurrentperturbationsproducenophaseresettingandwherethePRCflipsitssign.Sincethereareonlytwoneutralpoints,allPRCsofplanarneuraloscillatorsarebimodal.ThedegreeofbimodalityofaPRC,thatis,theratiobetweentheamplitudesofthedelayandadvancelobesofaPRC,canbesmoothlyadjustedwhenthebifurcationscenarioleadingtostableoscillatorybehaviorcombinesasaddlenodeofinvariantcircle(SNIC)andanAndronov-Hopfbifurcation(HB).


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 2 ISRN Computational Biology 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, 1 ( ) = 1 / − 1, where is the intrinsic period of oscillation of isolated neuron and 1 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).
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 [11][12][13]. 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 [15][16][17][18]. 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 [19][20][21].
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) [22][23][24]. 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 0.175 [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 [22-24, 26, 27]. Recent theoretical studies suggested that a SNIC bifurcation can also lead to bimodal PRCs [28][29][30].
In this paper we shown that (1) all PRCs are bimodal and (2) 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 (1) introducing the concept of neutral points, that is, phases at which external perturbations produce no phase resetting, and (2) 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.

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 F = ( 1 ( 1 ( ), 2 ( )), 2 ( 1 ( ), 2 ( ))) , where x = ( 1 , 2 ) are the fixed reference frame variables, which in the case of Morris-Lecar model [27,33] are the membrane potential, 1 , and the potassium conductance, 2 , 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 framex = ( , ) of the figurative point traveling along the LC, , than in the fixed frame x = ( 1 , 2 ) 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,x = x, is not unique [34][35][36]. However, the unit vector of the tangent direction is uniquely defined by k t = ( 1 , 2 ) /√ 2 1 + 2 2 and then we chose one of the two possible orientations of the normal vector to be k n = ( 2 , − 1 ) /√ 2 1 + 2 2 . In the fixed reference frame, x, the temporal evolution of small perturbations, x, from is described by the linearized equations of the dynamical system; that is, ( x)/ ≈ x, where = ( / )| is the Jacobian matrix evaluated along . Using the unitary transform U, the same perturbation with respect to the moving reference frame,x, is given by [29,34] 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 = 0, 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 x, 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, Δ = (Δ, 0) with respect to the fixed reference frame, x (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, = 0, 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, Figure 2: A stable limit cycle (LC) generates a closed trajectory in the fixed reference frame ( 1 , 2 ) . (a) A membrane potential perturbation Δ moves the figurative point 1 from the unperturbed LC to 1 . The tangent projection of the phase space perturbation Δ is (0) and represents an instantaneous geometric phase jump that corresponds to a temporal phase resetting. The normal projection of Δ is (0) 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 1 would have traveled to 2 along the unperturbed LC whereas the perturbed figurative point 1 reached 2 . 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 1 ( ). (b) In the absence of any coupling between the normal and tangent directions of the LC, that is, = 0, the only phase resetting is determined by . When the perturbation Δ is perpendicular to the LC, it results in that (0) = 0 and there is no phase resetting (points " " and " ") at these two neutral points.

ISRN Computational Biology
which has the analytic solution: 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.

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 ( 0 ) in (5) is zero. Since the normal perturbation does not contribute to the tangent displacement (because = 0), 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 = 0, then the maximum and the minimum of the membrane potential are the neutral points of the limit cycle.
If ̸ = 0, 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, ( ) = ( 0 ), is exactly balanced by the contribution of the normal relaxation processes through ( ) = ( 0 ). 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 ( = 0), 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 1 because the tangent projection of the perturbation, ( 0 ), 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 1 and one slow variable 2 ) 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 ( = 0), then the ISRN Computational Biology 5 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.

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 1 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 0.005 (in absolute units, which is ≈2.45% of the bias current SNIC at the SNIC bifurcation point) and a duration of 0.5% . 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, ≈ 2.45% SNIC stimulus amplitude has a very large effect (over 25% maximum phase resetting) close to the SNIC bifurcation point (Figure 3(A1)), whereas for a stimulus amplitude that is only about 4% 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 ( Ca = 1.0) 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 LP 2 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 LP 3 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.

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: (1) on the hyperpolarized stable steady state branch there is indeed a SNIC bifurcation point for low values of the bias current, and (2) 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 (LP 2 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 iPRC 1 ( ) and iPRC 2 ( ); that is, PRC( ) = 1 * (1− cos( ))+ 2 * sin( + ), where 1 and 2 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, 2 ≈ 0. 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, 1 ≈ 0. Since the amplitudes  ) there is a subcritical Andronov-Hopf and double-limit cycle bifurcations with stable oscillations emerging at 0 = 0.2419 (see the point LP 3 ) 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 0 = 0.2037 (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 10 −4 (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 (iPRC 1 ( )) and Type II PRC (iPRC 2 ( )). The ratio of the coefficient of iPRC 1 ( ) against iPRC 2 ( ) 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.
8 ISRN Computational Biology current perturbation (amplitude and duration), we computed their ratio 1 / 2 , which is insensitive to such variations (see Figure 4(d)). Our numerical simulations show that indeed the contribution of iPRC 1 ( ) is one order of magnitude larger than that of iPRC 2 ( ) 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 1 / 2 ≈ 5.7, we can actually use Figure 4(d) to estimate the range of the bias currents for which such an arbitrary classification is valid.

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 ( Ca = 1). 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 Ca . For example, for large values of calcium conductance, Ca , 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 LP 2 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 LP 1 and LP 2 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 peakto-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 1 = 0.1 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 Ca , require a relatively modest depolarizing bias current compensation to sustain oscillations. This is the reason why the bias current range is narrower for larger Ca 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 voltagegated 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)).

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 [41][42][43]. Therefore, the conductances of voltagegated 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  Figure 5: Bifurcation diagrams of ML model neuron for variable bias current at high (a) and low (b) calcium conductances, respectively. For large Ca , stable limit cycle oscillations emerge at low bias currents through a SNIC bifurcation (point LP 2 ) 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 LP 1 and LP 2 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.
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 [46][47][48], 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,[54][55][56] and references therein). Neural networks oscillate over more than three orders of magnitude and include delta (0.5-3 Hz), theta (3)(4)(5)(6)(7)(8), gamma , 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 [64][65][66]. 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.
The two nullclines of (A.1) are two-dimensional curves 1 ( 1 , 2 ) = 0 and 2 ( 1 , 2 ) = 0, 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.