Multiple and Complex Codimension-2 Bifurcations underlying Post-inhibitory Rebound Spike and Excitability Transition

Neuron exhibits nonlinear dynamics such as excitability transition and post-inhibitory rebound (PIR) spike related to bifurcations, which are associated with information processing, locomotor modulation, or brain disease. PIR spike is evoked by inhibitory stimulation instead of excitatory stimulation, which presents a challenge to the threshold concept. In the present paper, 7 codimension-2 or degenerate bifurcations related to 10 codimension-1 bifurcations are acquired in a neuronal model, which presents the bifurcations underlying the excitability transition and PIR spike. Type I excitability corresponds to saddle-node bifurcation on an invariant cycle (SNIC) bifurcation, and type II excitability to saddle-node (SN) bifurcation or sub-critical Hopf (SubH) bifurcation or sup-critical Hopf (SupH) bifurcation. The excitability transition from type I to II corresponds to the codimension-2 bifurcation, Saddle-Node Homoclinic orbit (SNHO) bifurcation, via which SNIC bifurcation terminates and meanwhile big homoclinic orbit (BHom) bifurcation and SN bifurcation emerge. A degenerate bifurcation via which BHom bifurcation terminates and fold limit cycle (LPC) bifurcation emerges is responsible for spiking transition from type I to II, and the roles of other codimension-2 bifurcations (Cusp, Bogdanov-Takens, and Bautin) are discussed. In addition, diﬀerent from the widely accepted viewpoint that PIR spike is mainly evoked near Hopf bifurcation rather than SNIC bifurcation, PIR spike is identiﬁed to be induced near SNIC or BHom or LPC bifurcations, and threshold curves resemble that of Hopf bifurcation. The complex bifurcations present comprehensive and deep understandings of excitability transition and PIR spike, which are helpful for the modulation to neural ﬁring activities and physiological functions.


Introduction
The bifurcations have played important roles in identifying the nonlinear dynamics of neuronal electrical activities, which are involved in the physiological functions of the neural system such as information processing and locomotor control and pathology or diseases of the brain [1][2][3][4][5][6][7][8][9]. The neuronal electrical activities include the steady state and firing, and the transition between steady state and firing or between different firing patterns can be induced by various modulations via different bifurcations [10][11][12][13][14][15][16][17][18][19][20]. For example, the steady state changes to firing via the sub-critical Hopf (SubH) bifurcation, or the sup-critical Hopf (SupH) bifurcation, or the saddle-node bifurcation on an invariant cycle (SNIC), or the saddle-node (SN) bifurcation [21,22]. Different bifurcations exhibit different dynamics when modulated by the external static stimulation [23][24][25][26]. For example, the firing frequency near SNIC point exhibits nearly zero value while near the SubH point manifests nearly fixed, nonzero value. Such different firing frequency responses have been used to in-terpret the type I and II excitabilities observed in the biological experiment [23]. For type I excitability, the steady state changes to firing with nearly-zero frequency, which then increases gradually with increasing stimulation strength. For type II excitability, the firing frequency emerges with nearly fixed, nonzero value. Therefore, type I and II excitabilities correspond to SNIC and SubH bifurcations related to equilibrium points, respectively. In biological experiment, type I and II excitabilities have been observed in many neurons, for example, the somatosensory neurons related to pathological pain, spinal sensory neurons, cortical pyramidal neurons, mesencephalic trigeminal neuronsgranule cells of dentate gyrus, sympathetic neurons, and hippocampal output neurons [27][28][29][30][31][32][33][34][35][36][37]. In addition, the type of spiking is determined by the firing frequency of spiking changed to steady state as the current strength decreases. Type I spiking exhibits nearly zero firing frequency around the disappearance of stable limit cycle, which corresponds to SNIC or big homoclinic orbit (BHom) bifurcation. Type II spiking exhibits fixed, non-zero firing frequency, which corresponds to the fold limit cycle (LPC) bifurcation or sup-critical Hopf (SupH) bifurcation. All these show that type I and II excitabilities/spikings are very important for both nonlinear dynamics and neuroscience.
Except for the firing frequency, type I and II excitabilities exhibit very different dynamics in multiple aspects [38,39]. For instance, the firing for type I and II excitabilities exhibit different phase response curves to the pulse stimulation [40]. Stimulated by noise, type I excitability exhibits weak coherence resonance while type II excitability manifests strong coherence resonance [26]. Network composed of neurons with type II excitability exhibits stronger synchronization than that composed of neurons with type I excitability [41,42]. Another important example is the post-inhibitory rebound (PIR) spike evoked by the inhibitory stimulation from the steady state. The PIR spike has been widely observed in the neurons such as the neonatal rat motoneurons [43], lamprey neurons [44], interneurons in mollusk clione [45], lateral pyloric neurons of the stomatogastric ganglion [46], and superior paraolivary nucleus in the auditory brainstem [47], which is involved in information processing and locomotor modulation. It is widely accepted that the PIR spike can be evoked from the steady state near Hopf bifurcation rather than near SNIC bifurcation [22]. The PIR spike presents a challenge to the concept of the threshold. In common viewpoint, the spike or action potential is evoked by the positive or excitatory stimulation which can induce the membrane potential to increase to run across the threshold, while a negative or inhibitory s-timulation cannot evoke a spike from the steady state. The PIR spike has been observed in the biological experiments on many neurons and related to special ionic currents such as the hyperpolarization active caution current (I h ) [48][49][50][51][52], which plays important roles in the physiological functions information processing and locomotor control of the neural system or brains disease such as depression [3]. In phase plane of a twodimensional neuronal model such as the famous Morris-Lecar (ML) model, which is widely used to identify the dynamics of type I and II excitabilities [53][54][55][56], the threshold curve across which a spike is evoked locates down-left, below, and right to the stable focus near the Hopf bifurcation [22]. Here the phase plane is plotted by membrane potential on the horizontal axis and gating variable on the vertical axis. The down left part of the threshold curve exhibits a negative slope. When stimulated by the inhibitory pulse, the phase trajectory of the response can run across the down-left part of the threshold curve firstly and then runs to right and increases to form a PIR spike. However, such geometrical characteristics of the threshold curve vanish for the case of SNIC in the ML model (for details please refer to the Appendix and Fig. 16 of body content in the present paper). The threshold curve, which is related to the stable manifold of the saddle, mainly locates right to the stable node near the SNIC and does not exhibit negative slope. Therefore, the down-left part of the threshold curve is responsible for the generation of PIR spike for Hopf bifurcation corresponding to type II excitability. However, no PIR spike is evoked due to the absence of down-left part of the threshold curve for the SNIC with type I excitability. Recently, I h current induces SNIC changed to Hopf bifurcation, which further implies that the PIR spike induced by I h current is related to the Hopf bifurcation [57]. More recently, in a three-dimensional neuronal model with I h current, the PIR spike appears near SNIC because the threshold surface changes for stronger I h current [58]. However, up to now, no PIR spike is identified to appear near SNIC in a two-dimensional model. In addition, the PIR spike was conjectured to appear near a big homoclinic orbit (BHom) bifurcation [22]. Unfortunately, no theoretical model and examples have been provided up to now. Therefore, an important question arises. Can the PIR spike appear near SNIC in a two-dimensional neuronal or near BHom bifurcation? The answer is very important for the nonlinear dynamics and neuroscience related to the PIR spike.
Due to both types of excitability are very important, the dynamics or bifurcations related to type I excitability and type II excitability, and the dynamics or bifurcations for the transition between excitability types have attracted much attention in both theoretical studies [59][60][61][62][63][64][65] and biological experiments [27][28][29][30][31][32][33][34][35][36][37]. In the biological experiments, the switch between type I excitability and type II excitability was identified and induced by multiple modulations such as the M-type potassium current, A-type current, optogenetic stimulation, and in vitro or in vivo-like conditions, and so on [28-30, 35, 57, 58, 64]. In theoretical models, the transition between type I excitability and type II excitability has been simulated in multiple theoretical models such as the ML model, I Na,p + I K model [23], Leech model, and so on, and related to the autapse, I h current, temperature, and so on. Furthermore, the bifurcations in two-parameter space [66][67][68][69][70] have been used to identify dynamical mechanism of excitability transition. The codimension-2 bifurcation point, Bogdanov-Takens (BT) bifurcation, which is related to codimension-1 bifurcations such as the SN bifurcation, homoclinic orbit (Hom) bifurcation, and Hopf bifurcation, is suggested to be related to switch between type I and II excitabilities. In some Refs. [5,62,65], the codimension-2 bifurcation point, Saddle-Node Homoclinic orbit (SNHO) bifurcation which is related to SNIC, SN, and BHom, is thought to be related to the transition between type I and II excitabilities. Therefore, the roles of codimension-2 bifurcations such as BT or SNHO bifurcations in the transition between type I and II excitabilities should be further identified.
In the present paper, the codimension-1 and -2 bifurcations underlying the PIR spike and the switch between excitability types are investigated in the I Na,p + I K model, which is a two-dimensional neuronal model without I h current. Such a model has been used to study the dynamics of type I and II excitabilities [23]. Firstly, the PIR spike is identified to appear near the SNIC, BHom, and LPC bifurcations, which is related to not only the type II excitability but also the type I excitability. The PIR spike near the SNIC/type I excitability presents a novel result, and the PIR spike near the BHom bifurcation presents a demonstration of the speculation proposed in a previous study [22]. Such results present the comprehensive viewpoints of the codimension-1 bifurcations and excitability types underlying the PIR spike. Secondly, the codimension-1 bifurcations underlying the excitability types are acquired. Type I excitability corresponds to the SNIC bifurcation, and type II excitability to the SN bifurcation right to the BHom, the SubH right to the BHom, the SubH right to the LPC, and the SupH. Type I spiking corresponds to the SNIC and BHom (left to the SN and SubH), and Type II spiking corresponds to the LPC (left to SubH) and SupH. Last, the codimension-2 and degenerate bifurcations underlying the PIR spike and switch between excitability types are acquired. There are multiple codimension-1 bifurcations and 7 codimension-2 or degenerate bifurcations in the I Na,p + I K model. Such bifurcations are more complex than those of the ML model and many other high-dimensional models such as the HH, Chay, Leech model, and so on. The excitability transition from type I to II corresponds to the codimension-2 bifurcation, Saddle-Node Homoclinic orbit (SNHO) bifurcation, via which the SNIC terminates, and meanwhile the BHom and SN bifurcations emerge. A degenerate bifurcation same as P 1 in Ref. [54], via which the BHom terminates and LPC emerges, is responsible for the spiking transition from type I to II. The roles of other codimension-2 bifurcations such as the Cusp bifurcation (CP), Bogdanov-Takens bifurcation (BT), and Bautin bifurcation (GH) are discussed. The novel bifurcation and excitability mechanisms for the PIR phenomenon and the comprehensive bifurcation mechanisms for the excitability transition are helpful for the modulations to the neural firing activities and physiological functions.
The rest of the present paper are organized as follows. Section 2 presents the model and method. The results are provided in Section 3. Section 4 presents the conclusion and discussion.

Model
The I Na,p + I K model [23] has been used to investigate the dynamics of excitability, which consists of a persistent Na + current, a relatively slower persistent K + current, and a leak current, reading as follows: where V is the membrane potential with unit mV and n denotes the unitless gating variable for K + channel. The parameter C represents the capacitance of a neuron and I denotes a constant current.
note the ionic currents of Na + , K + , and leakage ion respectively. E Na , E K , E L are the corresponding reversal potential and g Na , g K , g L are the corresponding maximal conductances. The functions m ∞ and n ∞ are steady state of the gating variables of Na + and K + , respectively, and γ is the time factor, where m ∞ (V ) = In the present paper, I is the bifurcation parameter and V 1/2 n the control parameter, and the other parameter values are k n = 7 mV, V 1/2 m = −20 mV, k m = 15 mV, C = 1 µF/cm 2 , E L = −79.42 mV, E K = −90 mV, E Na = 60 mV, g L = 8 mS/cm 2 , g K = 10 mS/cm 2 , g Na = 20 mS/cm 2 . And the unit of current is µA/cm 2 .

Threshold curve and the PIR spike
The threshold is an important characteristic of the neuronal electrical activity especially for the steady state. In common viewpoint, when an excitatory or depolarization stimulation is applied to a steady state, the membrane potential increases to go beyond a threshold to form an action potential. If the stimulation is weak, the membrane potential increases, but fails to reach the threshold, thus no action potential but subthreshold membrane potential appears. This is the common viewpoint of the voltage threshold.
In the present paper, the voltage threshold for the action potential evoked from the steady state is studied. For the I Na,p + I K model, if an initial value (V, n) for the Eq.(1) induces an action potential, such an initial value is called suprathreshold and labeled as blank in the phase plane. The initial value (V, n), which cannot induce an action potential, is called subthreshold and labeled with yellow in the phase plane. The border between the yellow and blank areas in the phase plane is the threshold curve. In the present paper, if the maximal value of the membrane potential is greater than 0 mV, an action potential is thought to appear. Obviously, the phase point of the steady state locates within the yellow area. The phase trajectory from the steady state to an action potential, which is evoked from the steady state, runs across the threshold curve, because a part of the action potential such as the peak locates within the blank area. And the phase point of the steady state and the subthreshold membrane potentials evoked from the steady state locate within the yellow area. Therefore, the threshold curve is very important to identify the generation of an action potential evoked from the steady state, especially for the PIR spike.
Except for the well-known excitatory stimulation, if an inhibitory stimulation can elicit an action potential or spike from the steady state, i.e. the trajectory of the electrical behavior following the stimulation run across the threshold curve, then the PIR spike is evoked. In the present paper, the relationship between the threshold curve and PIR spike is studied. It should be emphasized here that the phase plane is always plotted with the membrane potential V horizontally and gating variable n vertically, which is aimed to keep consistent with that in Ref. [22] and thus is convenient to make a comparison.

The inhibitory stimulation
In the present paper, we adopt inhibitory (i.e., negative) pulse stimulation, denoted by I stim , to the steady state to study the PIR spike. Adding I stim into the right hand of the first equation in Eq. (1), one has When initial value (V, n) is set to be the value of the steady state, Eq. (2) governs the response to the inhibitory pulse stimulation. Throughout the whole paper, I stim is a simple pulse with application time at 50 ms and termination time at 51.3 ms, i.e., duration fixed at 1.3 ms. And the pulse strength is labeled with A (< 0) µA/cm 2 .

Method
Fourth-order Runge-Kutta method is utilized to integrate the equations of the model with time step 0.001 ms. The bifurcation analysis is applied to Eq. (1) and performed by software XPPAUT [71], including oneparameter and two-parameter bifurcations.

The electrical behaviors
The I Na,p + I K model exhibits three different types of behaviors, i.e., the resting state, the period-1 firing, and the depolarization block. For example, at V 1/2 n = −29 mV, the representatives for the 3 behaviors are depicted in Fig. 1(a) (I = 3 µA/cm 2 , stable node corresponding to resting state), Fig. 1(b) (I = 100 µA/cm 2 , stable limit cycle corresponding to period-1 firing), and Fig. 1(c) (I = 240 µA/cm 2 , stable focus corresponding to depolarization) respectively. The membrane potential of the resting state and depolarization block are about −59.83 mV ( Fig. 1(a)) and about −18.98 mV ( Fig. 1(c)). In the present paper, the dynamics near the resting state and the period-1 firing and transition between them are studied.

Two-parameter bifurcation diagram
To acquire comprehensive view of bifurcation phenomena, bifurcations in the plane (I, V 1/2 n ) are acquired, as shown in Fig. 2. A groups of codimension-1 bifurcation curves locates in the left part and a sup-critical Hopf (SupH 2 , cyan) bifurcation locates at higher I value. The SupH 2 is related to the bifurcation from limit cycle to depolarization block, which is not the focus of the present paper. The codimension-1 bifurcation curves in left part are enlarged in Fig. 2(b), which is the focus of the present paper. Fig. 3(a), (b), (c), and (d) represent local enlargement of Fig. 2(b) around different codimension-2 or degenerate bifurcation points, where the gray region is for the mono-stable equilibrium, pink region for the stable limit cycle, blank region for the coexistence of stable equilibrium and stable limit cycle. There are 10 codimension-1 bifurcation curves, including SNIC 1 , B-Hom, SN 1 , SN 2 , SHom 1 , LPC, SubH, SupH 1 , SNIC 2 , and SHom 2 , respectively. Here SNIC 1 and SNIC 2 are for saddle-node bifurcation on an invariant circle, LPC for fold bifurcation of limit cycle, SubH for sub-critical Hopf, SupH 1 for sup-critical Hopf, SN 1 and SN 2 for saddle-node bifurcation, BHom and SHom 1 and SHom 2 for homoclinic bifurcation. There are 7 codimension-2 or degenerate bifurcation points, which are intersection points of some codimension-1 bifurcation curves. Four out of them are the well-known Saddle-Node Homoclinic orbit bifurcation (SNHO 1 , red solid circle), Bogdanov-Takens bifurcation (BT, black solid circle), Cusp bifurcation (CP, red solid circle), and Bautin bifurcation (GH, cyan solid circle), as shown in Fig. 2 Other three points locate within the big dashed circle (at the black arrow) shown in Fig. 2(b), which are two degenerate bifurcation points P 1 (Fig. 3(c)) and N  Fig. 2(b). The degenerate bifurcation point P 1 (red hollow star) is responsible for the change from BHom to LPC bifurcation; (d) In the dashed rectangle (at red arrow) of Fig. 3(c). The curve SHom 2 (purple dashed) connects points SNHO 2 (black hollow star) and a degenerate bifurcation point N (orange hollow star). The gray, pink, and blank region denotes the mono-stable equilibrium, stable limit cycle, and the coexistence.
Codimension-2 bifurcation point SNHO 1 Shown in Fig. 3(a) is the enlargement of Fig. 2(b) around the point SNHO 1 (red solid circle), which corresponds to the intersection of the codimension-1 bifurcation curves SN 1 (red solid line), BHom (blue dotted line), and SNIC 1 (green solid line). Via V 1/2 n decreasing across SNHO 1 , the curve SNIC 1 (green solid line) changes to curve SN 1 , and meanwhile the curve BHom emerges. Upper to the point SNHO 1 , the curve SNIC 1 separates exactly the steady state (gray) and spiking (pink) and there is no coexistence behavior. Lower to the SNHO 1 , the coexistence behavior (blank region) appears and locates between curves BHom and SN 1 . Therefore, SNHO 1 corresponds to the appearance of the BHom bifurcation and the coexistence behavior with decreasing V 1/2 n .
Codimension-2 bifurcation point BT As V 1/2 n further decreases, the SN 1 curve related to the SNHO 1 bifurcation point ( Fig. 3(a)) runs across the BT point and then the SubH curve (black dashed line) and SHom 1 curve (cyan dashed line) emerge from the BT point. Therefore, the BT bifurcation is responsible for the e-mergence of the SubH (black dashed line)and the SHom 1 (cyan dashed line) bifurcations, as shown in Fig. 3(b). The coexistence behavior locates between curves BHom and SN 1 upper to BT, whereas between curves BHom and SubH lower to BT. That is to say, the left boundary of the coexistence behavior (blank region) corresponds to curve BHom and does not change, whereas the right boundary changes from curve SN 1 to curve SubH.
Degenerate bifurcation point P 1 Fig. 3(c) shows the details in the big dashed circle depicted in Fig.  2(b). With decreasing V 1/2 n , the curve BHom (blue dotted) terminates and the curve LPC (magenta solid) appears via point P 1 (red star). Meanwhile, the curve SN 2 (dashed red) terminates and a saddle-node bifurcation on an invariant cycle (labeled as SNIC 2 , green dashed) curve appears through the P 1 point. Such a point P 1 has been reported in Morris-Lecar model ( Fig.  6 of [54]), which is thought to be a degenerate bifurcation point therein. Moreover, one can find that the left boundary of the coexistence behavior (blank region) changes from curve BHom to curve LPC via point P 1 , whereas the right boundary is still curve SubH.
Codimension-2 bifurcation point SNHO 2 and degenerate bifurcation point N Shown in Fig. 3(d) are the details of Fig. 3(c) within the dashed rectangle (corresponding to red arrow). With decreasing V 1/2 n , the SNIC 2 curve (green dashed line) changes to the SN 2 curve (red dashed) via the point SNHO 2 (black hollow star), meanwhile a SHom 2 curve (purple dashed line) appears. The curve SHom 2 (purple dashed) emerging from the SNHO 2 bifurcation and the SHom 1 curve (cyan dashed line) emerging from the point BT ( Fig.  3(b)) contact with each other to form a degenerate bifurcation, labeled with N (orange hollow star) in the present paper, as shown in Fig. 3(d). Such a point N corresponds to the disappearance (or appearance) of curves SHom 2 and SHom 1 . The dynamics around the SNHO 2 point and N point resemble those around the degenerate bifurcation point P 2 in Fig. 6 of Ref [54].
Codimension-2 bifurcation point GH With decreasing V 1/2 n , the SubH curve (black dashed line) emerging from the BT point (black solid circle, Fig. 2(b) and Fig. 3(b)) and the LPC curve (magenta solid line) emerging from the degenerate bifurcation point P 1 (red star, Fig. 3(c)) contact with each other at point GH (cyan solid circle) and then disappear, and meanwhile a sup-critical Hopf bifurcation curve (SupH 1 , black solid curve) appears, as shown in Fig. 2(b). The coexistence behavior locates between curves LPC (magenta solid) and SubH (black dashed) when V 1/2 n is upper to GH, whereas disappears when V 1/2 n is lower to GH. That is to say, GH is responsible for the disappearance of the coexistence behavior.
Codimension-2 bifurcation point CP The Cusp bifurcation (CP) is the intersection of codimension-1 bifurcation curves SN 1 and SN 2 , as shown in Fig. 2(b) and the insert figure therein, which corresponds to the disappearance of both SN 1 and SN 2 with decreasing V 1/2 n .
These bifurcation phenomena can be further well interpreted by one-parameter bifurcation diagram with respect to I at different V 1/2 n values, which are introduced in following subsection.
3.3 One-parameter bifurcation diagram and types of (spiking, excitability) In this subsection, we explain the bifurcation dynamics around above 7 points in two-parameter bifurcation plane by showing one-parameter bifurcation diagrams with respect to I. Meanwhile, the corresponding change of firing frequency is acquired to describe the spiking type and excitability type. In the present paper, to show the detailed excitability type as current I is increased and spiking type as current I is decreased, the notation (A, B) is adopted, where A is either I or II for spiking type, and B is either I or II for excitability type. The spiking type and excitability type here are consistent with those of Refs. [22,23]. Such a notation (A,B) can be conveniently used to describe the types of (spiking, excitability). For example, (I, II) means type I spiking and type II excitability.
3.3.1 One-parameter bifurcation diagrams around SNHO 1 : the excitability type transition When V 1/2 n = −29 (upper to SNHO 1 point), the S-NIC bifurcation with respect to I appears, as shown in Fig. 4(a1), and Fig. 4(a2) and (a3) represent the details of Fig. 4(a1) around SNIC 1 point (green solid circle). The equilibrium is depicted by the blue S-shaped curve, where the upper branch (short for UB) represents the unstable focus (long dashed line) and stable focus (solid blue line) separated by the sup-critical Hopf bifurcation point SupH 2 (black solid circle), the middle branch (MB) denotes the saddle (dotted line), and the lower branch (LB) is formed by the stable focus (thick blue) and stable node (thin blue). The extreme values of stable limit cycle bifurcated from the SupH 2 are denoted by the lower and upper red curves. The limit cycle contacts with the intersection point (green solid circle) of the MB and LB at I ≈ 3.03631 µA/cm 2 to form the SNIC 1 bifurcation (indicated by the vertical dashed line), which corresponds to the codimension-1 bifurcation curve SNIC 1 (green) in Fig. 2 and Fig. 3(a). The intersection of the UB and MB is for a saddle-node (SN) bifurcation point, denoted by SN 2 (red hollow circle), which corresponds to the codimension-1 bifurcation curve SN 2 (red dashed) in Fig. 2. The symbol FN (at red arrow) on the LB denotes the transition between the stable focus and stable node. The stable equilibrium on LB and the stable limit cycle correspond to the gray area and pink area in Fig. 3(a) respectively. No coexistence of the firing and stable node appears near the SNIC, which is consistent with the result shown in Fig.  2 or Fig. 3(a). Here the bifurcation point SupH 2 (black solid circle) on UB corresponds to the codimension-1 bifurcation curve SupH 2 (cyan) in Fig. 2(a).
When V 1/2 n = −29.8, which is lower to SNHO 1 point, the big homoclinic orbit (BHom) bifurcation emerges, labeled as BHom (at the black arrow) in Fig.  4(b1-b3), where Fig. 4(b2) and (b3) represents the details of Fig. 4(b1) around the red solid circle (labeled as SN 1 ). The BHom bifurcation occurs at I = I BHom (≈ 3.5204736 µA/cm 2 ), as indicated by the vertical dashed line with arrow, and is formed by the intersection between the limit cycle (red) and saddle on the M-B, and meanwhile the stable node on LB is within such a homoclinic loop (more details are shown in Fig. 5). The symbol SN 1 is the intersection between the LB and MB of the S-shaped curve, and denotes the saddle-node  Fig. 2. The equilibrium is depicted by the blue S-shaped curve, where the upper branch (short for UB) represents the unstable equilibrium (long dashed line) and stable focus (thick blue), the middle branch (MB) denotes the saddle (dotted line), and the lower branch (LB) is formed by the stable focus (thick blue) and stable node (thin blue). The intersection of the UB and MB is for a saddle-node (SN) bifurcation point, denoted by SN 2 (red hollow circle) in the present paper. The symbol FN (at red arrow) on the LB denotes the transition between the stable focus and stable node. The upper and lower red solid lines are for the extreme values of the stable limit cycle, which emerges from the sup-critical Hopf ( SupH 2 , black solid circle ) bifurcation. (a1) and its enlargements (a2-a3) correspond to V 1/2 n = −29 mV. The green solid circle, i.e., the intersection of the MB and LB at I ≈ 3.03631 µA/cm 2 and meanwhile the emergence of a limit cycle form the SNIC bifurcation (labeled as SNIC 1 ), as indicated by the vertical dashed line. (b1) and its enlargements (b2-b3) correspond to V 1/2 n = −29.8 mV. The limit cycle (red solid line) appears via the BHom bifurcation (dashed vertical line) at I BHom ≈ 3.5204736 µA/cm 2 , and the intersection between the LB and MB forms a saddle-node (SN 1 , red solid circle) bifurcation.
bifurcation and occurs at I ≈ 3.52159 µA/cm 2 . Here the bifurcation points BHom and SN 1 correspond to the codimension-1 bifurcation curves BHom (blue dotted) and SN 1 (red solid) in Fig. 2 or Fig. 3(a). One finds clearly that the BHom appears left to the SN 1 , and thus the coexistence of the firing and stable node appears between the BHom and SN 1 , as shown in Fig. 4(b2-b3), which is consistent with the coexistence region (blank) lower to the SNHO 1 in Fig. 3(a). The BHom bifurcation is right to FN point, i.e., the transition point between the stable focus (bold blue curve) and the stable node (thin blue curve). That is to say, the steady state left to and close to the BHom bifurcation corresponds to the stable node. Other symbols and plots are the same as those shown in Fig. 4(a1-a3) such as SupH 2 , which is not provided in the following figures due to such a bifurcation is not the focus of the present paper.
The results of Fig. 4(a3) and Fig. 4(b3) show that the SNHO 1 is responsible for the bifurcation transition from SNIC to BHom and the emergence of SN, which are consistent with the result shown in Fig. 2 or 3(a). Fig. 5 describes the dynamical behavior of the big homoclinic orbit (black) in the phase plane, which corresponds to the BHom bifurcation (at black arrow, I BHom ≈ 3.5204736 µA/cm 2 ) depicted in Fig. 4(b1-b3). Here the red line and blue line are for the nullclines dV /dt = 0 and dn/dt = 0, respectively, as depicted in Fig. 5(a). The intersection between both nullclines form three equilibrium points, a stable node (orange point), a saddle (half-filled point), and an unstable focus (hollow point), which correspond to the LB, MB, and UB in Fig. 4(b1-b3), respectively. The position relationship between the stable node and saddle is clearly shown in the insert figure of Fig. 5(a) and Fig. 5(c). Superposing the phase  Fig. 4(b1-b3)). (a)The V -nullcline dV /dt = 0 (red) and n-nullcline dn/dt = 0 (blue). Their intersection points are stable node (solid circle), saddle (half-filled circle), and unstable focus (hollow circle); (b) The big homoclinic orbit (BHom, black) with the nullclines; (c) The enlargement of (b) around the stable node and saddle point; (d) The further enlargement of (b) around saddle. The solid and dashed green line represent the stable and unstable manifold of the saddle, respectively. The black arrow indicates the direction for both the trajectory and manifolds of saddle.
trajectory of the BHom orbit onto Fig. 5(a) forms Fig.  5(b) and the details around the saddle are enlarged in Fig. 5(c). The black circle, i.e., the BHom orbit, goes through the saddle (half-filled point) to form a big loop surrounding the stable node (orange point) with anticlockwise (black arrow). The minimal membrane potential of the BHom orbit is lower than that of the stable node. Furthermore, the details of the homoclinic orbit near the saddle are enlarged again in Fig. 5(d), where the solid (dashed) green line represents the stable (unstable) manifold of the saddle. After starting from the saddle, the trajectory of the orbit runs along the unstable manifold (dashed green) of the saddle to form the action potential. Before returning to the saddle, the trajectory of the BHom orbit runs along the stable manifold (solid green line) of the saddle. However, for a common homoclinic orbit, or called small homoclinic (SHom) orbit, the node is outside of the SHom orbit, i.e. the minimal membrane potential of the SHom orbit is higher than that of the node. Fig. 6(a) is the change of firing frequency near the SNIC 1 bifurcation in Fig. 4(a2-a3), where the transition between the stable equilibrium (steady state) and stable limit cycle (spiking) occurs at the SNIC 1 .

Shown in
With increasing I, the stable steady state changes to spiking via SNIC 1 bifurcation and with nearly zero frequency, as shown by the red right arrow in Fig. 6(a). With decreasing I, the spiking changes to stable steady state also via SNIC 1 bifurcation and with nearly zero frequency, as shown by the down arrow in Fig. 6(a). Therefore, both the spiking and excitability belong to type I at V 1/2 n = −29. The types of (spiking, excitability) are labeled as (I, I). Shown in Fig. 6(b) is the change of the firing frequency near the BHom bifurcation depicted in Fig. 4(b2-b3). It can be found that the transition from spiking to steady state is via the BHom bifurcation and exhibits nearly zero frequency, as shown by the red arrow in Fig. 6(b). However, the transition from resting state to spiking occurs at the SN 1 bifurcation and manifests nonzero frequency, as shown by the green arrow in Fig. 6(b). Then it follows that the types of (spiking, excitability) are (II, I) at V 1/2 n = −29.8. Therefore, with decreasing V 1/2 n from upper side to lower side of SNHO 1 , the spiking remains unchanged and belongs to type I, whereas the excitability changes from type I to II.  Fig. 4(a2-a3)); (b) V 1/2 n = −29.8 mV (corresponding to Fig. 4(b2-b3)). Fig. 7: The one-parameter bifurcation diagram to explain the dynamics around the BT point in Fig. 2(b) or Fig.  3(b). The BHom bifurcation occurs at the black arrow (dashed vertical line). The LB exhibits stable focus (bold) around the BHom bifurcation point. (a1) and its enlargement (a2) correspond to V 1/2 n = −30.5 mV. (b1) and its enlargement (b2) correspond to V 1/2 n = −32.5 mV. A sub-critical Hopf (SubH, black hollow point) bifurcation appears on the LB. The unstable limit cycle (cyan dashed line) contacts with the saddle on the MB to form another homoclinic orbit bifurcation (SHom 1 , cyan point). Other plots and symbols are same with those in Fig. 4. Fig. 7(a1) corresponds to the one-parameter bifurcation diagrams at V 1/2 n = −30.5 mV (upper to BT) and its enlargement around SN 1 is further depicted in Fig.  7(a2), where the bifurcation structure and symbols are similar to those in Fig. 4(b2-b3). The only difference is that the BHom bifurcation point changes to locate left to the FN point, which means that it is stable focus (bold blue) close to the BHom bifurcation, as shown in Fig. 7(a2). The steady state left to the BHom becomes stable focus instead of stable node in Fig. 4(b3). As V 1/2 n further decreases to −32.5 mV (lower to BT), the BHom, SN 1 , and SN 2 still exist, as depicted in Fig.  7(b1). The obvious change is that the equilibrium point on LB exhibits a sub-critical Hopf (SubH, black hollow point) bifurcation and the unstable limit cycle (cyan dashed line) contacts with the saddle on the MB to form another homoclinic orbit bifurcation (labeled as SHom 1 , cyan point), as shown in Fig. 7(b2), which represents the details of Fig. 7(b1) around SN 1 . Here the bifurcation points SubH and SHom 1 correspond to the codimension-1 bifurcation curves SubH (black dashed) and SHom 1 (cyan dashed) in Fig. 3(b) respectively. Moreover, one sees clearly the coexistence for stable limit cycle and stable equilibrium locates between B-Hom and SN 1 in Fig. 7(a1-a2), whereas between BHom and SubH in Fig. 7(b1-b2). These results are consistent with those observed in Fig. 3(b).

One-parameter bifurcation diagrams around BT
Shown in Fig. 8 is the change of firing frequency at V 1/2 n = −32.5 mV (corresponding to Fig. 7(b1-b2)). One clearly sees that the transition from spiking to resting occurs at BHom with nearly zero firing frequency, as indicated by the red arrow. Therefore the type of spiking still corresponds to BHom and belongs to type I. Moreover, as indicated by the green arrow, it is via SubH bifurcation with a nearly fixed nonzero value that the resting changes to spiking, which thus shows type II excitability. Therefore, at V 1/2 n = −32.5 (lower to BT), the types of (spiking, excitability) are (I, II). Note that the types of (spiking, excitability) are also (I, II) in Fig. 6(b) (V 1/2 n = −30.5, upper to BT). Then one has that neither the spiking nor the excitability changes with V 1/2 n decreasing via BT, although the BT bifurcation induces the bifurcation related to equilibrium point changed from the SN 1 to the SubH.

One-parameter bifurcation diagrams around P 1 : the spiking type transition
The dynamics around degenerate point P 1 (Fig. 3(c)) can be explained by the one-parameter bifurcation diagrams shown in Fig. 9, where Fig. 9(a) corresponds to V 1/2 n = −33.2 mV (upper to point P 1 ), Fig. 9(b) to V 1/2 n = −33.30452 mV (corresponding to point P 1 ), and Fig. 9(c) to V 1/2 n = −33.3 mV (lower to point P 1 ). The LPC (magenta solid circle) bifurcation in Fig. 9(c) denotes the fold bifurcation of stable limit cycle, which means that, as I decreases, the stable limit cycle (red curves) and unstable limit cycle (orange dashed curves) coalesce and then disappear via LPC. The insert figure Fig. 9: The one-parameter bifurcation with respect to I at different V 1/2 n values around the degenerate bifurcation point P 1 in Fig. 3(c). Shown in dotted blue is the MB for the saddle and in dashed blue is the UB for the unstable focus. The orange (cyan) dashed line is for the unstable limit cycle. (a) V 1/2 n = −33.2 mV. The points SHom 1 (cyan) and SHom 2 (orange) are for the saddle homoclinic orbit bifurcation; (b) V 1/2 n = −33.2845 mV. The point P 1 (red hollow star) corresponds to that in Fig. 3(c); (c) V 1/2 n = −33.3 mV. The LPC (magenta solid circle) denotes the fold bifurcation of limit cycle and is formed by the intersection of stable limit cycle (red solid curves) and unstable limit cycle (orange dashed curves). The SNIC 2 (green hollow circle) denotes a SNIC bifurcation and is formed by that the disappearance of unstable limit cycle (orange dashed curves) occurs at the intersection of UB and MB. Other plots and symbols are the same as those in Figs. 4 and 7.
in Fig. 9(c) shows that the unstable limit cycle (orange dashed line) disappears at the intersection point of UB and MB, which thus denotes a SNIC bifurcation (labeled as SNIC 2 ). Other plots in Fig. 9 are the same as those in Fig. 7(b1-b2). Here the bifurcation points LPC and SNIC 2 correspond to the codimension-1 bifurcation curves LPC (magenta) and SNIC 2 (green dashed) in Fig. 3(c). The stable limit cycle (red curves) in Fig.  9(a) terminates via the BHom bifurcation, in Fig. 9(b) via the red star (the intersection of UB and MB), in Fig. 9(c) via the LPC bifurcation (magenta solid circle). Therefore, the red star in Fig. 9(b) corresponds to the point P 1 in Fig. 3(c) and represents the transition from BHom ( Fig. 9(a)) to LPC (Fig. 9(c)) with decreasing V 1/2 n . Meanwhile, via such point P 1 , SN 2 in Fig. 9(a) changes to SNIC 2 in Fig. 9(c). Besides, the coexistence behavior locates between BHom and SubH in Fig. 9(a) and between LPC and SubH in Fig. 9(c). These results are consistent with those observed in Fig.  3(c).  . The type of spiking now is determined by the firing frequency near the LPC bifurcation, which exhibits a nearly fixed, nonzero value. That is to say, the spiking belongs to type II. Besides, the excitability is still type II corresponding to the SubH bifurcation. Therefore, for V 1/2 n = −33.3 mV (lower to point P 1 ), types of (spiking, excitability) are (II,II). Note that the types of (spiking, excitability) at V 1/2 n = −32.5 (upper to point P 1 ) are (I,II), as depicted in Fig. 8. That is to say, as V 1/2 n decreases via the point P 1 , the spiking changes from type I to type II, meanwhile the excitability type keeps unchanged and still belongs to II. The results show that the degenerate bifurcation point P 1 plays an important role in modulating spiking between type I and type II.

One-parameter bifurcation diagrams around points N and SNHO 2
To further show the dynamics near the degenerate bifurcation point N in Fig. 3(d), one-parameter bifurcation diagrams with respect to I at different V 1/2 n values around the point N are shown in Fig. 11, which represents the local view of the bifurcation structure around the intersection of the MB and UB of equilibrium (blue) curve. When V 1/2 n = −33.3043 mV (upper to N), the SHom 2 (orange point) locates left to the SHom 1 (cyan point), as shown in Fig. 11(a). Here SHom 2 denotes a homoclinic orbit bifurcation, which is formed by the unstable limit cycle (orange dashed line) contacting with the saddle on the MB, and corresponds to the codimension-1 bifurcation curve SHom 2 (purple dashed) in Fig. 3(d). For V 1/2 n = −33.30452 mV (at point N), the SHom 2 and the SHom 1 contact with each other to form degenerate bifurcation N (orange hollow star), as shown in Fig. 11(b). When V 1/2 n = −33.3048 mV (lower to point N), both SHom 2 and SHom 1 disappear, as shown in Fig. 11(c). Therefore, the degenerate point N describes how the SHom 2 and SHom 1 coalesce and disappear, as depicted in Fig. 3(d). Since SHom 2 and SHom 1 are related with unstable limit cycle and unstable equilibrium, point N has no direct influences on the spiking or excitability transition.
Moreover, one can also explains the dynamics of SNHO 2 in Fig. 3(d) from the one-parameter bifurcation diagrams shown in Fig. 9(c) (V 1/2 n = −33.3, upper to SNHO 2 ) and Fig. 11(a) (V 1/2 n = −33.3043, lower to SNHO 2 ). The bifurcation structures in Fig. 9(c) and Fig. 11(a) are almost identical except for two pieces of difference. One striking difference is that the unstable limit cycle (orange dashed) terminates via the SNIC 2 bifurcation in former whereas via SHom 2 bifurcation in latter. The other one is that, the intersection of UB and MB denotes SNIC bifurcation (labeled as SNIC 2 ) in Fig. 9(c), whereas SN bifurcation (labeled as SN 2 ) in Fig. 11(a). Different from the SNHO 1 depicted in Fig. 3(a), the dynamical behaviors shown in Fig. 3(d) related to the 3 bifurcation curves through the SNHO 2 are unstable. More precisely to say, the curve SN 2 (red dashed) is related to the unstable node, the curve SNIC 2 (green dashed) related to the unstable node and unstable limit cycle, as depicted in Fig. 9(c) and Fig. 11(a), respectively, and the curve SHom 2 (purple dashed) related to the unstable limit cycle. Therefore, the SNHO 2 has no direct influence on the spiking or excitability transitions.

One-parameter bifurcation diagrams around GH
The one-parameter bifurcation diagrams shown in Fig.  12 can detail the dynamics of Bautin bifurcation (GH, cyan solid circle) in Fig. 2(b). When V 1/2 n = −34 mV (upper to GH), it exhibits LPC (magenta solid circle) and SubH (black hollow circle) bifurcation, and the coexistence behavior locates between LPC and Sub-H, as depicted in Fig. 12(a). As V 1/2 n decreases to be lower than GH, for example, at −40 mV, it exhibits sup-critical Hopf (SupH 1 , black solid circle) bifurcation, and no coexisting behavior appears, as depicted in Fig.  12(b). Here the bifurcation point SupH 1 corresponds to the codimension-1 bifurcation curve SupH 1 (black solid) in Fig. 2(b). Shown in Fig. 13 is the change of firing frequency at V 1/2 n = −40 (corresponding to Fig.  12(b)). The transitions between the steady state and spiking occur both via the SupH 1 bifurcation, where the firing frequency is a nearly fixed, non-zero value. Therefore, the types of (spiking, excitability) are still (II, II) at V 1/2 n = −40 (lower to GH). The GH point does not induce the transition between excitability type or spiking type when compared with those shown in Fig.  10 (V 1/2 n = −33.3, upper to GH). Fig. 12: One-parameter bifurcation diagrams to detail the dynamics of GH in Fig. 2(b). The equilibrium (blue) curve exhibits a monotonic characteristic. (a)V 1/2 n = −34 mV. The stable limit cycle (red solid) and the unstable limit cycle (dashed orange) contacts to form the LPC bifurcation (magenta solid circle). SubH (black hollow circle) denotes the sub-critical Hopf bifurcation; (b) V 1/2 n = −40 mV. SupH 1 (black solid circle) denotes the sup-critical Hopf bifurcation.

One-parameter bifurcation diagrams around CP
When V 1/2 n is upper to the CP, the equilibrium point curve exhibits a S-shaped curve and two SN bifurcation points SN 1 and SN 2 , as depicted in Figs. 4, 7, and 9. When V 1/2 n is lower to the CP, the equilibrium point curve becomes a monotonous characteristic and the SN bifurcation of the equilibrium points disappears, as depicted in Fig. 12. Therefore, the Cusp bifurcation (CP) determines the shape of the equilibrium point curve in (I, V ) bifurcation plane, and also corresponds to the disappearance of the SN bifurcation. The point CP has nothing to do with the spiking or excitability transitions due to that it is unrelated to the bifurcation of limit cycle.

The PIR spike for SNIC
The negative stimulation can induce spike from the steady state near the SNIC bifurcation depicted in Fig.  4(a1-a3), whereas no spike is evoked near the SNIC for the ML model in Ref. [22]. Therefore, our result presents a novel viewpoint of the post-inhibitory rebound (PIR) spike. Shown in Fig. 14  The threshold curve for the steady state near the SNIC bifurcation The PIR spike or the subthreshold membrane potentials can be explained by the equilibrium points and the threshold dynamics in the phase plane (V, n). Shown in Fig. 15(a) are the nullclines dV /dt = 0 (red) and dn/dt = 0 (blue) for V 1/2 n = −29 mV and I = 3.03 µA/cm 2 (left and close to I SNIC ≈ 3.03631 µA/cm 2 in Fig. 4(a2-a3)). Their intersection points denote the stable node (orange solid point), the saddle (half-filled point), and the unstable focus (hollow point). The stable node is a global stable attractor and corresponds to the resting state in Fig. 14. The insert figure shows the details around the stable node and the saddle.
For Eq. (1), the initial value of (V, n) to evoke a spike locates within the blank area, and the initial value to evoke subthreshold membrane potentials within the Fig. 14: The response of the steady state to an inhibitory pulse current for V 1/2 n = −29 mV and I = 3.03 µA/cm 2 (left and close to the SNIC bifurcation in Fig.  4(a2-a3)). The inhibitory pulse current with strength A = −115.56 µA/cm 2 (lower panel, black) can evoke a PIR spike (upper panel, black), while A = −115.55 µA/cm 2 (lower panel, red dotted) can induce not a PIR spike but the subthreshold membrane potential (upper panel, red dotted).
yellow area, as shown in Fig. 15(b1-b3) and Fig. 15(c1-c3). The border between the yellow and blank areas is the threshold curve. The square denotes the phase point at the termination time of the inhibitory pulse stimulation.
The trajectory for A = −115.55 µA/cm 2 , which corresponds to the subthreshold membrane potential (red dotted) in Fig. 14 and is evoked from the stable node (orange solid point), wholly locates within the yellow area, as shown by the red dotted curve in Fig.  15(b1). Its details around the square and saddle are further enlarged in Fig. 15(b2) and Fig. 15(b3) respectively. The trajectory at the termination time (square) of the inhibitory stimulation locates within the yellow area, which shows that the inhibitory stimulation with A = −115.55 µA/cm 2 cannot induce the trajectory (red dotted) to run across the threshold curve, as shown in Fig. 15(b2). And the trajectory near the saddle (halffilled point) locates left to the stable manifold (stable green line) of the saddle, i.e. the membrane potential lower than that of the saddle.
However, the trajectory for A = −115.56 µA/cm 2 , which corresponds to the PIR spike (black) in Fig. 14, exhibits very different dynamics, as shown by the black curve in Fig. 15(c1), (c2), and (c3), where the latter two represent the enlargement of the former one near the square and saddle respectively. The inhibitory stimulation with A = −115.56 µA/cm 2 induces the trajectory (black) to run across the down left part of the threshold curve with a negative slope, i.e. go from the yellow area to blank area along the direction from up-right to down-left, as shown in Fig. 15(c2). The trajectory near Fig. 15: The dynamics for the PIR spike or the subthreshold membrane potentials shown in Fig. 14. (a) The stable node (orange solid point), saddle (half-filled point), and unstable focus (hollow point) correspond to the intersection points between nullclines dV /dt = 0 (red) and dn/dt = 0 (blue). Inset shows the details around the node and saddle; (b1) The phase trajectory (red dotted) for A = −115.55 µA/cm 2 corresponding to the subthreshold membrane potential (red dotted) in Fig. 14; (b2) The enlargement of (b1) near the square (i.e., the phase point at the termination of the pulse stimulation); (b3) The enlargement of (b1) near the saddle; (c1) The phase trajectory (black) for A = −115.56 µA/cm 2 corresponding to the PIR spike (black) in Fig. 14; (c2) The enlargement of (c1) near the square; (c3) The enlargement of (c1) near saddle. The blank and yellow area represent the collection of the initial values to induce a spike and not a spike but threshold behaviors, respectively. The border between the yellow and blanket areas forms the threshold curve. The arrow indicates the flowing direction of the trajectories. The solid (dashed) green line in (b3) and (c3) denotes the stable (unstable) manifold of saddle. the saddle locates right to the stable manifold (green solid line) of the saddle and runs along the unstable manifold (green dashed line) of the saddle, as shown in Fig. 15(c3). Therefore, the right part of unstable manifold takes the role to increase the membrane potential near the saddle. Then the membrane potential further increases to form the PIR spike, as shown in Fig. 15(c1). To sum up, it is the part of the threshold curve downleft to the stable node, which exhibits a negative slope, that is responsible for the PIR spike evoked from the stable node near the SNIC. Once the inhibitory stimulation drives the trajectory to run across such part of the threshold curve at the termination time of the stimulation, the PIR spike is evoked, otherwise not. Besides, the part of the threshold curve near the saddle coincides exactly with the stable manifold (solid green) of the saddle, which shows that the dynamics of saddle plays important roles in the generation of the PIR spike. The shape of the threshold curve for the SNIC in the present paper resembles that of Hopf bifurcation in Ref. [22] to a large extent.

Different from the traditional viewpoint that no PIR spike for SNIC bifurcation in ML model
To further show the difference of the threshold curve for the SNIC between the I Na,p + I K model in the present paper and the well-known ML model, some results related to Ref [22] are acquired, as shown in Fig. 16 and Fig. A1 of Appendix. The stable node, the saddle, and the unstable focus are labeled with the orange solid point, half-filled point, and hollow point, respectively, and the threshold curve is the boundary between the yellow and blank areas, as shown in Fig. 16(a). The details near the saddle (solid and dashed curve represen- Fig. 17: The n-nullclines (solid line) exhibit different slope in phase plane (V, n) for I Na,p +I K model (red)and ML model (black). The solid circle denotes stable node. The parameters are same with those in Fig. 15 and Fig.  16 respectively. t the stable and unstable manifold of the saddle) are further depicted in Fig. 16(b), which shows that the threshold curve coincides with the stable manifold of the saddle. The whole trajectory evoked by A = −200 µA/cm 2 (red dotted) and A = −500 µA/cm 2 (black) locate in the yellow area, as shown in Fig. 16(c). The square denotes the phase point at the termination time of the stimulation. Even the stimulation is very strong (A = −500 µA/cm 2 ), no PIR spike but subthreshold membrane potential (black) is evoked. Such a result has been widely accepted. The equations, parameter values, bifurcation diagram, and phase portrait of the ML model, please see the Appendix.
Except for the different shapes of threshold curve between I Na,p + I K model and ML model, it can be found that the slope (dn/dV ) of n-nullcine around the stable node (solid circle) exhibits different slopes, as shown in Fig. 17. The slope of n-nullcine for I Na,p + I K  Fig. 4(b1-b3). A = −3.79 µA/cm 2 (black) and A = −3.78 µA/cm 2 (red dotted) denote the two critical strengths to and not to evoke the PIR spike; (b) The resting state (stable focus) is near BHom bifurcation shown in Fig. 7(b1-b2). The two critical strengths are A = −3.84 µA/cm 2 (black) and A = −3.83 µA/cm 2 (red dotted). Fig. 19: The dynamics of the PIR spike for the BHom bifurcation shown in Fig. 18. The left panel describes the phase trajectory of the PIR spike in Fig. 18(a) and the right panel corresponds to that in Fig. 18(b). The solid orange point denotes the stable equilibrium, half-filled point for the saddle, and hollow point for the unstable focus. The hollow square represents the phase point at the termination time of the stimulation. The yellow area denotes for the initial values to induce subthreshold behaviors, whereas the blank area for the initial values to induce spike behaviors. Their boundary forms the threshold curve for the spike. (a1)-(b1): The global view of the phase trajectory. Inset in (a1) shows the position relationship of node and saddle; (a2)-(b2): The enlargement of (a1)-(b1) near the squares, respectively; (a3)-(b3): The enlargement of (a1)-(b1) near the saddle, respectively, where the solid (dashed) green line represents the stable (unstable) manifold of the saddle. model (red) is larger than that of ML model (black), which influences the vector fields around the n-nullcine. The value of slope of n-nullcine is a possible important factor to determine the PIR spike or not for a same type of bifurcation, which should be further studied in future. Fig. 18 shows the PIR spike evoked from the steady state near the BHom bifurcation. For the stable node at I = 3.52 µA/cm 2 for V 1/2 n = −29.8 mV, which is left and close to the BHom (I BHom ≈ 3.5204736 µA/cm 2 ) shown in Fig. 4(b1-b3), the inhibitory stimulation with strength A = −3.79 µA/cm 2 (lower black) can induce PIR spike (upper black), as shown in Fig. 18(a), while the inhibitory stimulation with A = −3.78 µA/cm 2 (lower red dotted) induces not PIR spike but subthreshold membrane potentials (upper red dotted).

The PIR spike near BHom
For the stable focus at I = 5.75 µA/cm 2 (V 1/2 n = −32.5 mV), which is left and close to the BHom (I BHom ≈ 5.75239 µA/cm 2 ) shown in Fig. 7(b1-b2), the inhibitory stimulation with strength A = −3.84 µA/cm 2 (lower black) can induce PIR spike (upper black), as shown in Fig. 18(b), while the inhibitory stimulation with A = −3.83 µA/cm 2 (lower red dotted) induces the subthreshold membrane potentials (upper red dotted). The membrane potential finally recovers to the steady state. The main difference between Fig. 18(a) and Fig. 18(b) is the absence of damping oscillations for the former and the presence for the latter, as depicted by the black arrows, which is due to that the steady state is stable node for the former and stable focus for the latter.
The phase trajectory for A = −3.79 µA/cm 2 is illustrated in left panel of Fig. 19, which correspond to the PIR spike in Fig. 18(a) and is evoked from the stable node. The inset in Fig. 19(a1) shows clearly the position relationship between node (orange point) and saddle (half filled point). Shown in right panel of Fig. 19 is the phase trajectory for A = −3.84 µA/cm 2 , which correspond to the PIR spike in Fig. 18(b) and is evoked from stable focus (orange point). As depicted in Fig. 19, these two trajectories for PIR spikes related to BHom exhibit dynamics similar to that of the SNIC (Fig. 15). The down-left part of both threshold curves exhibit a negative slope, as shown in Fig. 19(a1) and (b1). The enlargement of the threshold curve with the negative slope near the square (the termination time of the inhibitory stimulation) is shown in Fig. 19(a2) and (b2). The inhibitory stimulation can induce the trajectory run across the threshold curve from up-right to down-left. The enlargement of the trajectory of the PIR spike (black) near the saddle (half-filled point) is shown in Fig. 19(a3) and (b3). The trajectory locates right to the stable manifold (solid green curve) and unstable manifold (dashed green curve), resembling that of the SNIC. The hollow circle in Fig. 19(a1) and Fig.  19(b1) represents the unstable focus. Fig. 20 is the PIR spike near the LPC bifurcation (I LPC ≈ 6.64876 µA/cm 2 ) in Fig. 9(c). Fixed I = 6.64 µA/cm 2 (left and close to LPC bifurcation), the system has a unique equilibrium, i.e., a stable focus (orange point). The inhibitory stimulation with strength A = −4.1 µA/cm 2 (black, lower panel) can evoke a PIR spike (black, upper panel) and A = −4.09 µA/cm 2 (red dotted, lower panel) induce not PIR spike but subthreshold membrane potential (red dotted, upper panel), as shown in Fig. 20(a). The phase trajectory (black) of the PIR spike exhibits anti-clockwise (black arrow), as depicted in Fig. 20(b). The details around the focus (orange point) and the threshold curve at the square (the termination of the stimulation) are shown in Fig.  20(c) and (d), respectively. The part of the threshold curve locating down-left to the stable focus (orange point) exhibits a negative slope, and the inhibitory stimulation can induce the trajectory run across the threshold curve from up-right to down-left. The square represents the termination time of the stimulation, and damping oscillations appear before recovering to the steady state.

Shown in
3.5 Summary of the bifurcations underlying the excitability/spiking transition and the PIR spike In summary, the type I and type II excitabilities are separated by the horizontal line across the codimension-2 bifurcation SNHO 1 (red solid point), and type I and type II spiking are separated by the horizontal line across the degenerate bifurcation point P 1 (white hollow star), as shown in Fig. 21. Therefore, the plane (I, V 1/2 n ) is divided into 3 regions. From top to bottom, the types of (spiking, excitability) for the three regions are (I, I), (I, II), and (II, II), which are called case-1, -2, and -3 respectively in the present paper. Furthermore, either of case-2 and case-3 contains two subcases. The two subcases for case-2 are separated by the thin line across the codimension-2 bifurcation BT (black solid point), and are labeled as case-2 1 and case-2 2. And the two subcases for case-3 are separated by the thin line across the codimension-2 bifurcation GH (cyan solid point), and are labeled as case-3 1 and case-3 2.
The bifurcations related to types of (spiking, excitability) are (SNIC, SNIC), (BHom, SN 1 ), (BHom,  The relationship between the codimension-2 (or degenerate) bifurcation points and the types of spiking/excitability in two-parameter plane (I,V 1/2 n ). According to the distributions of spiking/excitability type and the related bifurcations, the whole plane is divided into five cases by four dashed horizontal lines, which are labeled with case-1, -2 1, -2 2, -3 1, -3 2, respectively, from top to bottom. According to the neuronal dynamical behaviors, the whole plane is divided into three regions, which are marked with gray (mono-stable equilibrium), pink (stable limit cycle), and blank (coexistence of stable equilibrium and stable limit cycle), respectively. The other plots are same with those shown in Figs. 2-3. , which lies between the steady state (gray) and spiking (pink), appears between the B-Hom and SN 1 bifurcation for the case-2 1, between the BHom and SubH bifurcations for case-2 2, and between the LPC and SubH bifurcations for case-3 1, as shown in Fig. 21 and Table 1.
For the steady state near SNIC, BHom, and LPC bifurcation points, the PIR spike is identified to be evoked. The PIR spike near the SupH bifurcation is not studied in the present paper due to the small amplitude of the firing, as shown in Fig. 12(b).

Conclusion and Discussion
Bifurcation has been successfully used to identify dynamics and physiological roles of the neuronal electrical activities. In the present paper, the multiple and complex codimension-1 and codimension-2 bifurcations underlying the PIR spike induced by the inhibitory stimulation and the transition between type I and type II excitabilities are acquired, which present novel results and comprehensive viewpoints in the following aspects.
Firstly, a novel example for the PIR spike, which appears near the SNIC bifurcation corresponding to type I excitability, is provided, and the speculation that the PIR spike can appear near the BHom bifurcation in Ref [22] is verified in a two-dimensional neuronal model. The PIR spike appearing near the SNIC bifurcation is a novel result, and the PIR spike near the BHom bifurcation presents a demonstration to the previous spec-ulation. Combined the PIR spike appearing near the LPC bifurcation corresponding to type II excitability and near the Hopf bifurcation in the previous Refs [22,23,[54][55][56], the comprehensive bifurcations for the PIR spike are acquired. Different from the previous viewpoint that the PIR spike appears near the bifurcation corresponding to type II excitability instead of type I excitability [22], the PIR spike is identified to appear near the SNIC bifurcation related to type I excitability in the present paper, which extends the emerging condition for the PIR spike.
Secondly, a novel threshold curve across which the PIR spike is evoked from the stable node is presented for the SNIC or BHom bifurcations in the I Na,p + I K model, which is different from the threshold curve of the ML model. When the membrane potential is adopted as horizontal ordinate of the phase plane, the striking geometrical feature is that a part of the threshold curve locates down-left to the stable node and exhibits a negative slope, which resembles the threshold curve for the Hopf bifurcation (a part of the threshold curve locates down-left to the stable focus). However, such a part of threshold curve cannot be found for the SNIC bifurcation in the famous ML model. Last, the comprehensive viewpoint of the codimension-2 (or degenerate) bifurcations underlying excitability and spiking transition is presented. Although the I Na,p + I K model is very simple, for example, only two variables and two-dimensional, seven codimension-2 or degenerate bifurcations and 11 codimension-1 bifurcation curves are acquired. Type I excitability corresponds to the SNIC bifurcation, while type II excitability corresponds to the SN (right to BHom), SubH (right to B-Hom or LPC), and SupH bifurcations. Therefore, type I excitability changes to type II excitability via the codimension-2 bifurcation, SNHO, wherein the SNIC bifurcation terminates and meanwhile the SN and B-Hom bifurcations emerge. Coexistence of the steady s-tate and firing appears between the BHom and SN bifurcations, which leads to that the firing frequency at the SN point is non-zero. Therefore, the SN bifurcation corresponds to type II excitability. Type I spiking corresponds to the SNIC and BHom bifurcation points, while type II spiking corresponds to the LPC and SupH bifurcations. Therefore, a degenerate bifurcation point P 1 , via which the BHom bifurcation terminates and the LPC bifurcation emerges, is responsible for the spiking transition from type I to type II. In addition, the I Na,p +I K model presents bifurcations in two-parameter space more complex than those in the ML model, H-H model, and many other neuronal model with highdimension and multiple variables [13,[54][55][56][57][58][66][67][68][69][70].
In future, the excitability transition and codimension-2 or degenerate bifurcations induced by other parameters or modulations for the I Na,p + I K model, the excitability transition and codimension-2 or degenerate bifurcations in other neuronal model (such as degenerate points P 2 and c 1 in Ref. [54]), and the cause for the distinct threshold curves for the PIR spike evoked near the SNIC bifurcation between the ML model and I Na,p + I K model should be further studied, for example, the influences of slope of n-nullcline on the different shapes of threshold curve. Through these studies, the deep understanding and comprehensive viewpoints for the excitability transition and PIR spike and the related bifurcations can be acquired, which are helpful for the identification of the nonlinear dynamics and physiological roles of the PIR spike and excitability transition in neuroscience. Availability of data and material Data will be made available on reasonable request.

Declarrations
Code availability Not applicable.
Ethics approval Not applicable.
Consent to participate Not applicable.

Consent for publication Not applicable.
The Morris-Lecar (ML) model in Refs. [3,53] reads as follows where V denotes the membrane potential with unit mV and n represents the unitless gating variable of K + ionic channel. C is for the capacitance of a neuron and I for a constant current. I K = −g K n(V −E K ), I Ca −g Ca m∞(V )(V −E Ca ), and I L = −g L (V − E L ) denote the ionic currents of K + , Ca 2+ , leakage ions, respectively. The parameters E K , E Ca , and E L are the corresponding reversal potentials, and the parameters g K , g Ca , and g L are corresponding the maximal conductances. The functions m∞ and n∞ are the steady state of the gating variables Ca 2+ and K + , respectively, γ is the time factor, where The parameter values are as follows: V 1 = −1.2 mV, V 2 = 18 mV, V 3 = 12 mV, V 4 = 17.4 mV, E L = −60 mV, E K = −84 mV, E Ca = 120 mV, g L = 2 µS/cm 2 , g K = 8 µS/cm 2 , g Ca = 4 µS/cm 2 , and C = 20 µF/cm 2 . The parameter I is considered as bifurcation parameter.
The bifurcation diagram for the SNIC bifurcation (I SNIC ≈ 39.96 µA/cm 2 ) is shown in Fig. A.1(a), and the equilibrium points and nullclines at I =39.95 µA/cm 2 (left and close to SNIC) in phase plane are shown Fig. A.1(b). The phase portrait at I = 39.95 µA/cm 2 . The intersection points between the nullclines dV /dt = 0 (red) and dn/dt = 0 (blue) denote the equilibrium points, i.e., the stable node (orange point), saddle (half-filled point), and unstable focus (hollow point), which are the same as those shown in Fig. 16 in body text.