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

Authors: Sorinel A. Oprisan
Journal: ISRN Computational Biology
 (2013)

---

## Abstract

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


---

## Body

## 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 f is continuously tunable by adjusting a bias current I and such neurons can fire with arbitrarily low frequencies, or Class 2, that is, discontinuous f-I 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 periodPi, 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, F1(φ)=P1/Pi-1, where Pi is the intrinsic period of oscillation of isolated neuron and P1 is the transiently modified duration of the current cycle due to a perturbation applied at the stimulus time ts or phase φ=ts/Pi (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).The stable membrane potential oscillatory solution of Morris-Lecar model neuron (continuous line in panel (a)) has an intrinsic periodPi, which can be transiently changed to P1 (dashed line in panel (a)) due to a brief (inhibitory) perturbation applied at the stimulus time ts (vertical arrow) or stimulus phase φ=ts/Pi. Equally, spaces phases, 5%Pi apart, are marked by solid circles along the unperturbed membrane potential trace (continuous line). The first order phase resetting curve (b), F1(φ)=P1/Pi-1, represents the normalized change in the transient period of oscillation and is always bimodal. The amplitude, M, of the delay (F1>0) lobe is very large compared to the amplitude, m, of the advance (F1<0) lobe when the ML oscillator is very close to its SNIC bifurcation point (see inset in panel (b)).
(a)
(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 [11–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–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–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–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 Figure1(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–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.


## 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 areF=(f1(x1(t),x2(t)),f2(x1(t),x2(t)))T, where x=(x1,x2)T are the fixed reference frame variables, which in the case of Morris-Lecar model [27, 33] are the membrane potential, x1, and the potassium conductance, x2, and the superscript T 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~=(xt,xn)T of the figurative point traveling along the LC, γ, than in the fixed frame x=(x1,x2)T of the original system of differential equations (see Figure 2(a)). The choice of the unitary matrix U that switches between the reference frames, that is, x~=Ux, is not unique [34–36]. However, the unit vector of the tangent direction xt is uniquely defined by vt=(f1,f2)T/f12+f22 and then we chose one of the two possible orientations of the normal vector to be vn=(f2,-f1)T/f12+f22.A stable limit cycle (LC) generates a closed trajectory in the fixed reference frame(x1,x2)T. (a) A membrane potential perturbation Δ moves the figurative point A1 from the unperturbed LC to B1. The tangent projection of the phase space perturbation Δ is xt(0) and represents an instantaneous geometric phase jump that corresponds to a temporal phase resetting. The normal projection of Δ is xn(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 tϵ, the figurative point A1 would have traveled to A2 along the unperturbed LC whereas the perturbed figurative point B1 reached B2. 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 F1(φ). (b) In the absence of any coupling between the normal and tangent directions of the LC, that is, λtn=0, the only phase resetting is determined by λtt. When the perturbation Δ is perpendicular to the LC, it results in that xt(0)=0 and there is no phase resetting (points “a” and “c”) at these two neutral points.
(a)
(b)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, d(δx)/dt≈Aδx, where Aij=(∂fi/∂xj)|γ 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]
(1)d(δx~)dt=d(Uδx)dt=(U′U-1+UAU-1)δx~=1f12+f22(λttλtnλntλnn)δx~=Λδx~,
where λtt=(f12f1,x1+f22f2,x2+f1f2(f1,x2+f2,x1)), λtn=(f22-f12)(f1,x2+f2,x1) + 2f1f2(f1,x1-f2,x2), λnt=0, and λnn=f12f2,x2+f22f1,x1-f1f2(f1,x2+f2,x1) (we also corrected a typo in λtn of [34]). The coefficient λαβ measures the effect along the direction α for a perturbation applied along the direction β with respect to γ. For example, the coefficient λtn 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 λnt=0, which means that a perturbation tangent to γ never induces a normal displacement of the figurative point.While the Jacobian matrixA 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)T 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 U:
(2)Δ~=Δf12+f22(f1,f2)T.
Equation (2) gives the initial conditions with which we integrated the perturbation equations (1). Since for any planar dynamical systems, λnt=0, the differential equation for normal perturbations, that is, xn′=λnn(t)xn, has the analytic solution:
(3)xn(t,t0)=xn(t0)exp(∫t0t‍λnn(s)ds),
where the initial condition is determined by the injected current of amplitude Δ; that is, xn(t0)=(Δf2(t)/f1(t)2+f2(t)2)|t=t0. 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 λtt coefficient) and normal perturbations (through λtn); that is,
(4)xt′=λttxt+λtnxn=λttxt+λtnxn(t0)exp(∫t0t‍λnn(s)ds),
which has the analytic solution:
(5)xt(t)=exp(∫t0t‍λtt(s)ds)×((∫t0s‍(λnn(u)-λtt(u))du)xt(t0)+xn(t0)ssdddds×∫t0t‍(λtn(s)exp(∫t0s‍(λnn(u)-λtt(u))du)ds)),
where xt(t0)=(Δf1(t)/f1(t)2+f2(t)2)|t=t0.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.


## 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 areF=(f1(x1(t),x2(t)),f2(x1(t),x2(t)))T, where x=(x1,x2)T are the fixed reference frame variables, which in the case of Morris-Lecar model [27, 33] are the membrane potential, x1, and the potassium conductance, x2, and the superscript T 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~=(xt,xn)T of the figurative point traveling along the LC, γ, than in the fixed frame x=(x1,x2)T of the original system of differential equations (see Figure 2(a)). The choice of the unitary matrix U that switches between the reference frames, that is, x~=Ux, is not unique [34–36]. However, the unit vector of the tangent direction xt is uniquely defined by vt=(f1,f2)T/f12+f22 and then we chose one of the two possible orientations of the normal vector to be vn=(f2,-f1)T/f12+f22.A stable limit cycle (LC) generates a closed trajectory in the fixed reference frame(x1,x2)T. (a) A membrane potential perturbation Δ moves the figurative point A1 from the unperturbed LC to B1. The tangent projection of the phase space perturbation Δ is xt(0) and represents an instantaneous geometric phase jump that corresponds to a temporal phase resetting. The normal projection of Δ is xn(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 tϵ, the figurative point A1 would have traveled to A2 along the unperturbed LC whereas the perturbed figurative point B1 reached B2. 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 F1(φ). (b) In the absence of any coupling between the normal and tangent directions of the LC, that is, λtn=0, the only phase resetting is determined by λtt. When the perturbation Δ is perpendicular to the LC, it results in that xt(0)=0 and there is no phase resetting (points “a” and “c”) at these two neutral points.
(a)
(b)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, d(δx)/dt≈Aδx, where Aij=(∂fi/∂xj)|γ 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]
(1)d(δx~)dt=d(Uδx)dt=(U′U-1+UAU-1)δx~=1f12+f22(λttλtnλntλnn)δx~=Λδx~,
where λtt=(f12f1,x1+f22f2,x2+f1f2(f1,x2+f2,x1)), λtn=(f22-f12)(f1,x2+f2,x1) + 2f1f2(f1,x1-f2,x2), λnt=0, and λnn=f12f2,x2+f22f1,x1-f1f2(f1,x2+f2,x1) (we also corrected a typo in λtn of [34]). The coefficient λαβ measures the effect along the direction α for a perturbation applied along the direction β with respect to γ. For example, the coefficient λtn 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 λnt=0, which means that a perturbation tangent to γ never induces a normal displacement of the figurative point.While the Jacobian matrixA 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)T 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 U:
(2)Δ~=Δf12+f22(f1,f2)T.
Equation (2) gives the initial conditions with which we integrated the perturbation equations (1). Since for any planar dynamical systems, λnt=0, the differential equation for normal perturbations, that is, xn′=λnn(t)xn, has the analytic solution:
(3)xn(t,t0)=xn(t0)exp(∫t0t‍λnn(s)ds),
where the initial condition is determined by the injected current of amplitude Δ; that is, xn(t0)=(Δf2(t)/f1(t)2+f2(t)2)|t=t0. 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 λtt coefficient) and normal perturbations (through λtn); that is,
(4)xt′=λttxt+λtnxn=λttxt+λtnxn(t0)exp(∫t0t‍λnn(s)ds),
which has the analytic solution:
(5)xt(t)=exp(∫t0t‍λtt(s)ds)×((∫t0s‍(λnn(u)-λtt(u))du)xt(t0)+xn(t0)ssdddds×∫t0t‍(λtn(s)exp(∫t0s‍(λnn(u)-λtt(u))du)ds)),
where xt(t0)=(Δf1(t)/f1(t)2+f2(t)2)|t=t0.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 λtn vanishes. This means that the only contribution to the PRC is due to tangent perturbations through λtt. 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 xt(t0) in (5) is zero. Since the normal perturbation does not contribute to the tangent displacement (because λtn=0), the points where the brief (infinitesimal) voltage perturbation is normal to γ are the neutral points, that is, the points “a” and “c” in Figure 2(b). In other words, if λtn=0, then the maximum and the minimum of the membrane potential are the neutral points of the limit cycle.Ifλtn≠0, the exact positions (or phases) of the two neutral points a′ and b′ in Figure 2(b) can be determined from the requirement that the contribution of the pure tangential perturbations to (5), that is, xt′(t)=λttxt(t0), is exactly balanced by the contribution of the normal relaxation processes through xt′(t)=λtnxn(t0). 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 (λtn=0), then the only contribution to the PRC is due to λtt. On the other hand, λtt contribution vanishes at local minima or maxima of the membrane potential x1 because the tangent projection of the perturbation, xt(t0), 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 x1 and one slow variable x2) 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 (λtn=0), 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 a′ 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 durationP1 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 ISNIC at the SNIC bifurcation point) and a duration of 0.5%Pi.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)(I0=0.08334), 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, I0=0.09 (B1), I0=0.1 (C1), and I0=0.2 (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).
(a)
(b)
(c)
(d)Close to the SNIC bifurcation point (see Figure3(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, a≈2.45%ISNIC 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.Thef-I 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 (gCa=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 LP2 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 LP3 in Figure 4(b). To further support our conjecture that all PRCs are bimodal, we numerically computed the ratio m/M of the PRC lobes (see Figure 1(b)). Very close to the SNIC bifurcation point the amplitude m of the smaller lobe is a few orders of magnitude smaller than the magnitude M of the other lobe (Figure 4(c)), which justifies the name “unimodal” for Type I PRC, associated with SNIC bifurcations.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 conductancegCa=1 (b) shows a clear SNIC bifurcation with the fold characterized by the two turning points LP1 and LP2, 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 I0>0.15 in panel (b)) there is a subcritical Andronov-Hopf and double-limit cycle bifurcations with stable oscillations emerging at I0=0.2419 (see the point LP3) 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 I0=0.2037 (HP point in panel (b)). The ratio m/M 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 (iPRC1(φ)) and Type II PRC (iPRC2(φ)). The ratio of the coefficient of iPRC1(φ) against iPRC2(φ) 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.
(a)
(b)
(c)
(d)


### 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 Figure4(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 Figure4(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 (LP2 point in Figure 4(b)), which produces a PRC that can be approximated with the “unimodal” Type I expression
(6)iPRC1(φ)=1-cos(φ)
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
(7)iPRC2(φ)=sin(φ+α),
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 iPRC1(φ) and iPRC2(φ); that is, PRC(φ)=a1*(1-cos(φ))+a2*sin(φ+α), where a1 and a2 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, a2≈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,a1≈0. Since the amplitudes a1 and a2 depend on the characteristics of the injected current perturbation (amplitude and duration), we computed their ratio a1/a2, which is insensitive to such variations (see Figure 4(d)). Our numerical simulations show that indeed the contribution of iPRC1(φ) is one order of magnitude larger than that of iPRC2(φ) 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 a1/a2≈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.


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

The complete bifurcation diagram shown in Figure4(b) was obtained for fixed calcium conductance (gCa=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 gCa. For example, for large values of calcium conductance, gCa, 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 LP2 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 LP1 and LP2 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 V1=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 gCa, require a relatively modest depolarizing bias current compensation to sustain oscillations. This is the reason why the bias current range is narrower for larger gCa 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)).Bifurcation diagrams of ML model neuron for variable bias current at high (a) and low (b) calcium conductances, respectively. For largegCa, stable limit cycle oscillations emerge at low bias currents through a SNIC bifurcation (point LP2) 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 LP1 and LP2 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.
(a)
(b)


## 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 λtn vanishes. This means that the only contribution to the PRC is due to tangent perturbations through λtt. 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 xt(t0) in (5) is zero. Since the normal perturbation does not contribute to the tangent displacement (because λtn=0), the points where the brief (infinitesimal) voltage perturbation is normal to γ are the neutral points, that is, the points “a” and “c” in Figure 2(b). In other words, if λtn=0, then the maximum and the minimum of the membrane potential are the neutral points of the limit cycle.Ifλtn≠0, the exact positions (or phases) of the two neutral points a′ and b′ in Figure 2(b) can be determined from the requirement that the contribution of the pure tangential perturbations to (5), that is, xt′(t)=λttxt(t0), is exactly balanced by the contribution of the normal relaxation processes through xt′(t)=λtnxn(t0). 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 (λtn=0), then the only contribution to the PRC is due to λtt. On the other hand, λtt contribution vanishes at local minima or maxima of the membrane potential x1 because the tangent projection of the perturbation, xt(t0), 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 x1 and one slow variable x2) 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 (λtn=0), 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 a′ 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 durationP1 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 ISNIC at the SNIC bifurcation point) and a duration of 0.5%Pi.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)(I0=0.08334), 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, I0=0.09 (B1), I0=0.1 (C1), and I0=0.2 (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).
(a)
(b)
(c)
(d)Close to the SNIC bifurcation point (see Figure3(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, a≈2.45%ISNIC 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.Thef-I 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 (gCa=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 LP2 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 LP3 in Figure 4(b). To further support our conjecture that all PRCs are bimodal, we numerically computed the ratio m/M of the PRC lobes (see Figure 1(b)). Very close to the SNIC bifurcation point the amplitude m of the smaller lobe is a few orders of magnitude smaller than the magnitude M of the other lobe (Figure 4(c)), which justifies the name “unimodal” for Type I PRC, associated with SNIC bifurcations.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 conductancegCa=1 (b) shows a clear SNIC bifurcation with the fold characterized by the two turning points LP1 and LP2, 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 I0>0.15 in panel (b)) there is a subcritical Andronov-Hopf and double-limit cycle bifurcations with stable oscillations emerging at I0=0.2419 (see the point LP3) 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 I0=0.2037 (HP point in panel (b)). The ratio m/M 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 (iPRC1(φ)) and Type II PRC (iPRC2(φ)). The ratio of the coefficient of iPRC1(φ) against iPRC2(φ) 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.
(a)
(b)
(c)
(d)


## 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 Figure4(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 Figure4(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 (LP2 point in Figure 4(b)), which produces a PRC that can be approximated with the “unimodal” Type I expression
(6)iPRC1(φ)=1-cos(φ)
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
(7)iPRC2(φ)=sin(φ+α),
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 iPRC1(φ) and iPRC2(φ); that is, PRC(φ)=a1*(1-cos(φ))+a2*sin(φ+α), where a1 and a2 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, a2≈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,a1≈0. Since the amplitudes a1 and a2 depend on the characteristics of the injected current perturbation (amplitude and duration), we computed their ratio a1/a2, which is insensitive to such variations (see Figure 4(d)). Our numerical simulations show that indeed the contribution of iPRC1(φ) is one order of magnitude larger than that of iPRC2(φ) 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 a1/a2≈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.


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

The complete bifurcation diagram shown in Figure4(b) was obtained for fixed calcium conductance (gCa=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 gCa. For example, for large values of calcium conductance, gCa, 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 LP2 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 LP1 and LP2 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 V1=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 gCa, require a relatively modest depolarizing bias current compensation to sustain oscillations. This is the reason why the bias current range is narrower for larger gCa 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)).Bifurcation diagrams of ML model neuron for variable bias current at high (a) and low (b) calcium conductances, respectively. For largegCa, stable limit cycle oscillations emerge at low bias currents through a SNIC bifurcation (point LP2) 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 LP1 and LP2 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.
(a)
(b)


## 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 [41–43]. 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 Pi [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% Pi) [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 [46–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–56] 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 [64–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.



---

*Source: 230571-2013-12-12.xml*