Recent developments and future perspectives in aging and macrophage immunometabolism

1.   Introduction

Nonlinear dynamics has played very important roles in identifying dynamics of neural electrical activities [10,29,32], which are involved in the information processing, locomotion control, cognitive functions, and brain disease [8,19,31]. These activities mainly include the steady state such as the resting state and the firing behavior composed of action potentials or spikes, for example, the spiking and bursting, which has been widely investigated with the conceptions of bifurcations and chaos [6,16,18,23,33]. For example, by modulation to some physiological parameters, the steady state changes to firing behavior via the subcritical or supercritical Hopf bifurcation point, or the saddle-node bifurcation on an invariant cycle (SNIC), or the saddle-node bifurcation (SN) point, or the firing behavior changes to the steady state via the supercritical Hopf bifurcation point, or the fold bifurcation of limit cycle, or the SNIC, or the bifurcation of homoclinic (Hom) orbit [16,18]. The Hom bifurcation usually contains the bifurcation of big homoclinic orbit (BHom) and small homoclinic orbit(SHom, i.e., the common homoclinic orbit), which are named here according to the size of homoclinic loop. In addition, the SNIC, or SN, or Hom are built relationship to type Ⅰ excitability, which means that the firing generates with nearly zero frequency, and the Hopf bifurcation to type Ⅱ excitability, which corresponds to the firing with a nearly fixed frequency [16,18]. All these show that nonlinear dynamics is effective to identify the dynamics and physiological functions of the neural electronic behaviors.

The dynamics for an action potential or spike evoked by external stimulation from steady state has been an important issue in both neuroscience and nonlinear dynamics [15,16,17], which has been related to the nonlinear concept, threshold [15,16,17]. Well known that, a positive, or called excitatory, or called depolarization stimulation pulse with suprathreshold strength can induce membrane potential increased to reach a voltage threshold to elicit an action potential while with subthreshold strength can not evoke a spike but subthreshold membrane potential. In nonlinear dynamics, the threshold has been well explained with the trajectory or manifold in the phase space combined with the bifurcations as well as the types of excitability. For example, Fig. 14 in Ref. [16] describes how the action potential is evoked by excitatory stimulation. For a stable node near SNIC with type Ⅰ excitability [16,17,23], the threshold curve corresponds to the stable manifold of the saddle, which is right to the stable node and has a positive slope, i.e. the membrane potential of the saddle is larger than the stable node. The suprathreshold stimulation induces the membrane potential or phase trajectory increased to run across the stable manifold of the saddle to form an action potential, as depicted in Fig. 14 (left) in Ref. [16]. For the stable focus near Hopf bifurcation with type Ⅱ excitability [16,17,23], the threshold sets locate on three sides of stable focus, left, bottom, and right sides, as depicted in Fig. 14 (right) in Ref. [16]. Similarly to that of SNIC, it is also the part right to the stable focus that acts the threshold for an action potential induced by the excitatory stimulation. Generally speaking, for the excitatory stimulation, it is the right part of threshold curve or sets, which locates right to the steady state, are responsible for the generation of an action potential. In Ref. [28], the threshold curves for multiple two-dimensional models are acquired, however, the bifurcations for these models are not provided. The threshold curve for other bifurcations such as Hom awaits to be studied.

To the contrary, for the negative (or called inhibitory or hyperpolarization) stimulation, the left part of the threshold sets for the Hopf bifurcation is very important for the generation of an action potential, i.e. the post-inhibitory rebound (PIR) spike, which is depicted in Fig. 14 (right) of Ref. [16]. Therein the inhibitory stimulation induces the trajectory to run across the left part of the threshold sets, and then the trajectory locates outside of the threshold sets and rotates in anti-clockwise. After the trajectory approaches the right part of the threshold sets, the membrane potential increases to form an action potential. Therefore, it is the left part of threshold sets, locating left to the stable focus, are responsible for the generation of PIR spike. The PIR spike presents that inhibitory stimulation can enhance neural firing activity [2], which is different from the traditional viewpoint that the inhibitory modulations always suppress the electronic activity. It has been widely accepted that the PIR spike can be evoked from focus near the Hopf bifurcation point with type Ⅱ excitability instead of node near the SNIC with type Ⅰ excitability [16], due to the absence of the part of threshold curve left to the steady state for type Ⅰ excitability. In neurophysiology, the PIR spike has been build relationship to special ionic currents, mainly the hyperpolarization-active caution ($I_h$) current [1] or A-type potassium current [24]. The PIR spike identified in the nervous system with $I_h$ current is involved in the locomotor rhythms [3,26], "short-term memory" mechanism in the lateral pyloric neuron of the stomatogastric ganglion [11]. Recently, $I_h$ current is identified to induce the SN bifurcation changed to the Hopf bifurcation [35], which seems to follow that the PIR spike induced by $I_h$ current in neurophysiology and by Hopf bifurcation in nonlinear dynamics are consistent with each other. Therefore, it is often thought that the PIR spike in a neuron model is mainly related to the stable focus near the Hopf bifurcation with type Ⅱ excitability in nonlinear dynamics [2,16] or $I_h$ current in neurophysiology [1,3,11,26]. However, in Ref. [9], although the PIR spike (or burst, or spiking, or bursting) is simulated, their relationships to the bifurcations or threshold are not presented.

However, the condition for the PIR spike is extended to the SNIC in system with $I_h$ current in Ref. [12] and to canards of the fold saddle-type in Ref. [22]. In addition, in the system without $I_h$ current, the PIR spike is identified to appear near the SN bifurcation [18], which is build a relationship to a Bogdanov-Takens bifurcation, i.e., the intersection point between the SN, Hopf, and Hom bifurcation [16,18,25,34]. The results show that the non-Hopf bifurcation or non-$I_h$ current mechanism for the PIR has attracted increasing attention. Reviewing the literatures, one interesting conjecture for the PIR spike evoked from the steady state for the type Ⅰ excitability near a big homoclinic (BHom) bifurcation is depicted in Fig. 15 of Ref. [16]. The most striking characteristic is that the threshold curve for an action potential is around the stable steady state from the left, below, and right sides of the steady state, which resembles the threshold sets of the Hopf bifurcation with type Ⅱ excitability to a large extent. Unfortunately, no a theoretical model, or the detailed dynamics of the bifurcations and the PIR spikes are presented. Bifurcation of BHom is investigated in multiple theoretical models such as the Morris-Lecar model with autapse, modified FitzHugh-Nagumo (FHN) model, and $\rm {I_{Na}+I_{K}}$ model [16,18,34], and is related to a codimension-2 bifurcation, saddle-node Homoclinic orbit bifurcation (SNHO) (the intersection point of the codimension-1 bifurcation curve for the BHom, SNIC, and SN). Except for the type Ⅰ excitability, excitability transition, bifurcation itself, equilibrium points, and phase trajectory related to BHom bifurcation [16,18,34], the threshold curve and the PIR spike for the BHom bifurcation have not been studied.

In the present paper, the threshold curve and PIR spike for the BHom bifurcation are studied in a planar dynamical system without $I_h$ current. Firstly, the mono-stable node and stable node coexisting with spiking, which locates at the opposite sides of the BHom bifurcation point, are acquired. The PIR spike can be evoked from the stable node in both cases of the mono-stability and the coexistence, and the threshold curve surrounds the steady state from the left, bottom, and right sides. The left part of the threshold curve resembles that of the Hopf bifurcation but differs from that of the SNIC, and the right part resembles that of SNIC, i.e., coincides exactly with the stable manifold of the saddle, but differs from that of the Hopf bifurcation. Therefore, the threshold curve for the mono-stable node near the BHom bifurcation partially resembles but partially differs from those of the Hopf bifurcation and SNIC. The detailed dynamics of the trajectory of the action potential or spiking induced by both the inhibitory and excitatory stimulations can be well explained with the threshold curve, especially the intersection point between the trajectory and the threshold curve. Once the trajectory runs across the threshold curve during stimulation, the action potential is evoked, otherwise not. Secondly, the mono-stable foucs and the coexisted stable focus with spiking also locates on the opposite sides of the BHom bifurcation, and their threshold curves are qualitatively similar to those of the node, except for the damping oscillations following the action potential. Last, multiple codimension-1 bifurcation curves such as the BHom curve and one codimension-2 bifurcation point, the saddle-node Homoclinic orbit bifurcation point (SNHO), are acquired, which present the parameter region for the BHom bifurcation and the generation of the PIR spike. The results present the threshold curve for the BHom bifurcation, which is a novel case different from those of the well-known SNIC and Hopf bifurcation, and the PIR spike evoked from both focus and node near the BHom bifurcation, which extend the condition for the PIR spike from mainly Hopf bifurcation with type Ⅱ excitability or the $I_h$ current to the BHom bifurcation with type Ⅰ excitability unrelated to the $I_h$ current. Such results represent a comprehensive viewpoint to the threshold for the BHom bifurcation, which extends the concept of the threshold and condition for the PIR spike.

The rest of the paper is organized as follows. Section 2 shows the model and methods. The main results are displayed in Section 3. Section 4 gives the conclusions and discussions.

2.   Models and methods

2.1.   Model

The FitzHugh-Nagumo (FHN) model was created by reducing the Hodgkin-Huxley model to two dimensions. This reduced model successfully imitates the generation of action potential and therefore attracts a lot of attention, and has been studied by analytical, numerical and experimental methods [13,21,14]. In the present paper, we study the action potential induced by pulse current with a modified FHN model, which is Eq. (15) in Ref. [16] and reads as

where the variable $V$ mimics the membrane potential, the variable $w$ mimics activation of an outward current, and $u$ is treated as the bifurcation parameter. The parameter $\epsilon$ = 1 and the function $s(w) = \frac{b}{1+e^{(c-w)/d}}$ is an $S$-shaped function, where $b$ = 1.3, $d$ = 0.05 and $c$ is the control parameter. The model is dimensionless.

2.2.   Threshold curve

In the present paper, the threshold curve is acquired to study the action potential evoked from the steady state. The $(V, w)$ phase plane is divided into two parts by the threshold curve. One part denotes the collection of phase points whose values $(V, w)$ being as the initial values of the FHN model (i.e., Eq. (1)) can induce an action potential or spike. Such a part is labeled with blanket and is called suprathreshold area for convenience in the present paper. The other part, marked with yellow, represents the collection of phase points as the initial values of the FHN model to induce not action potential but the subthreshold potential, which is called subthreshold area. The border between these two parts is the threshold curve in $(V, w)$ plane. A suprathreshold stimulation can induce that the phase trajectory goes away from the steady state locating within the yellow region, runs across the threshold curve and into the blanket area, and at last forms a spike, which is the basic process to evoke a spike from the steady state. A subthreshold stimulation cannot induce the trajectory to run across the threshold curve from the subthreshold (yellow) area to the suprathreshold (blanket) area. In the present paper, action potential or spike refers to the maximal value of the membrane potential greater than 1.

2.3.   Stimulations to evoke spike or PIR spike

In the present paper, to study how the spike or PIR spike is evoked, we adopt a single pulse stimulation $I_{stim}$ and apply it to the steady state near the BHom bifurcation. That is to say, adding the pulse stimulation $I_{stim}$ into the right hand side of the first equation in Eq. (1) and meanwhile setting the initial values to be the steady state, then different stimulations induce different responses. The pulse duration is fixed to be 1, and the pulse strength is labeled with $A$. Since the spike is evoked from the stationary behavior of the steady state, the specific application time of $I_{stim}$ does not influence the generation of spike or not. For convenience, $I_{stim}$ always starts at time $t$ = 10 as representative in the present paper.

2.4.   Method

Fourth-order Runge-Kutta method is utilized to integrate the FHN model with time step $0.001$. The bifurcation is performed by software XPPAUT [7].

3.   Results

The results are divided into 3 parts. Firstly, a novel threshold curve for big homoclinic (BHom) bifurcation related to the stable node is acquired. The parameter value $c$ = $-$0.55 is chosen as representative. For the stable node left and close to BHom, the stimulation, regardless of excitatory and inhibitory one, can evoke single spike including PIR spike. When the stable node locates right and close to the BHom, spiking is evoked instead of single spike. The dynamical mechanisms for the evoked spike and spiking are obtained by comparing the phase trajectory and the threshold curve. Secondly, when $c$ = $-$0.4, the stable equilibrium near BHom bifurcation now is stable focus instead of stable node. Results similar to the stable node can be obtained. Finally, different bifurcations with respect to $u$ at different $c$ values and double-parameter bifurcations in ($u$, $c$) plane are acquired. Multiple codimension-1 bifurcation curves, such as BHom curve, and a codimension-2 bifurcation related to BHom bifurcation, the saddle-node Homoclinic orbit bifurcation point, are acquired. Different behaviors in the regions separated by bifurcations curves and PIR spike region around the BHom curve are acquired.

3.1.   Threshold for the stable node near BHom

3.1.1.   Bifurcations with respect to $u$ and the BHom bifurcation

The bifurcation diagram of equilibrium point (blue) and stable limit cycle (red) is shown in Fig. 1(a). Here the $S$-shaped blue line is formed by lower branch (LB), middle branch (MB), and upper branch (UB). The LB (thick solid blue line) denotes the stable node, MB (short dashed blue) represents the saddle, and UB is composed of unstable focus (long dashed blue line) and stable focus (thin solid blue line) corresponding to depolarization block. The symbol SN (black solid circle), i.e., the intersection point of LB and MB, denotes a saddle-node bifurcation at $u_{SN}$ $\approx$ $-$ 1.02. The symbol H (black hollow circle) represents a supercritical Hopf bifurcation via which a stable limit cycle emerges. Vmax (Vmin) denotes the maximal (minimal) value of the stable limit cycle. As $u$ decreases to $u_{BHom}$ $\approx$ $-$1.099400401984 (black arrow), a BHom bifurcation occurs. The stable limit cycle becomes a BHom orbit to the saddle on the MB. The insert shows the dynamics on UB at the red arrow, which are the stable node (blue solid line), the extreme values of stable limit cycle (green solid line), a saddle-node bifurcation point (red solid circle), a Hopf bifurcation point (red hollow circle) and a common Hom orbit bifurcation (half-filled circle, sometimes called small Hom orbit (SHom)). Due to these bifurcations are far from the BHom bifurcation, we do not pay attention to these bifurcations.

As $u$ decreased to $u_{BHom}$, the firing period (blue) increases almost to infinity and the firing frequency (red) decreases to nearly zero, as shown in Fig. 1(b), which shows that BHom bifurcation exhibits type Ⅰ excitability.

At the bifurcation point $u$ = $u_{BHom}$, the nullclines $dV/dt$ = 0 (red curve) and $dw/dt$ = 0 (blue curve), the BHom orbit (black loop with anti-clockwise arrow pointing), the stable node (red solid circle), the saddle (half-filled circle), and the unstable focus (red hollow circle) are shown in Fig. 2(a). This homoclinic orbit is so big that it can contain not only the unstable focus, but also the stable node. Therefore, the membrane potential of the left part of the orbit is lower than that of the stable node (red solid circle). Results similar to Fig. 2(a) please refer to Fig. 32 or 39 in Ref. [16]. The details of Fig. 2(a) near the saddle (half-filled circle) are shown in Fig. 2(b). The saddle (half-filled circle) locates on the orbit, and the solid and dashed green line represent the stable and unstable manifolds of the saddle, respectively. Therefore, the BHom orbit begins from the saddle, then runs to right and along the unstable manifold (dashed green), gets far from the saddle with increasing $V$ to form the action potential, and at last returns to the saddle exactly along the stable manifold (solid green) of the saddle from the bottom. Different from a BHom orbit which is big enough to contain a stable steady state, SHom orbit is relatively small and does not contain a stable steady state [4,16,18,20].

3.1.2.   Threshold for the mono-stable node (left and close to BHom)

For $u$ = $-$1.12 ($<u_{BHom}$), the system exhibits a mono-stable node on the LB, which corresponds to the resting state, as shown in Fig. 3(a). The membrane potential of stable node is labeled with $V_n$ and $V_n$ $\approx$ $-$1.00502342630403. The border between the yellow and blanket areas is the threshold curve, as shown in Fig. 3(b), where the blanket region is for the suprathreshold area and the yellow region for the subthreshold area. The solid red circle, half-filled circle, hollow red circle represent the stable node, saddle, and unstable focus, respectively. The membrane potential of saddle is labeled with $V_s$ and $V_s$ $\approx$ $-$0.703981477599643 for $u$ = $-$1.12($<u_{BHom}$). Such a shape of the threshold curve surrounds the stable node from the left, bottom, and right sides, which resembles that of Hopf bifurcation point. However, it differs from that of the SNIC due to that there does not exist the part of the threshold curve down-left to stable node near SNIC (Fig. 14 (left panel) in Ref. [16]). Therefore, the threshold curve for the stable node near the BHom bifurcation is a novel case different from that for SNIC.

A pulse stimulation with $A$ = 0.29 (red, lower panel) can not evoke a spike but with $A$ = 0.3 (black, lower panel) can induce a spike, as depicted in Fig. 4(a). For $A$ = 0.3, the membrane potential at the termination time (black arrow) of the stimulation is larger than $V_s$ (the purple dashed line). After the stimulation, the membrane potential further increases to form the spike. For $A$ = 0.29, the membrane potential at the termination time (black arrow) of the stimulation is lower than $V_s$. After the stimulation, the membrane potential decreases, therefore, no spike is induced. The membrane potential of saddle $V_s$ separates the suprathreshold (black) and subthreshold (red) potential, which suggests that the threshold for the spike is related to the saddle.

The generation of a spike (black curve) for $A$ = 0.3 (black) or subthreshold membrane potential (red curve) for $A$ = 0.29 (red) can be explained by the threshold curve, as shown in Fig. 4(b)-(d), wherein Fig. 4(c) and (d) are the enlargement of Fig. 4(b). The red trajectory for $A$ = 0.29 or subthreshold membrane potential locates within the subthreshold area (yellow), while the black trajectory containing an action potential runs across the threshold curve and locates within both subthreshold (yellow) and suprathreshold (blanket) areas. Such difference is induced by the different locations of the phase point (square) at the termination time of the pulse stimulation, as shown in Fig. 4(c). For $A$ = 0.29, the red hollow square is still within the yellow area, which shows that the phase trajectory beginning from the stable node (red solid circle) does not run across the threshold curve during the stimulation. However, for $A$ = 0.30, the black hollow square is within the blanket area, which shows that the phase trajectory beginning from the stable node (red solid circle) runs across the threshold curve during the stimulation. After the stimulation, the black trajectory approaches the right unstable manifold (dashed green line) of the saddle (half-filled circle) along the direction parallel to the stable manifold of the saddle (solid green line) firstly, then turns right and runs far away from the saddle along the right unstable manifold (dashed green line) of the saddle, as shown in Fig. 4(d), moves to right to a large extent to form the action potential, at last returns to the saddle, as shown in Fig. 4(b). After the stimulation, the red trajectory for $A$ = 0.29 approaches the unstable manifold left to the saddle first, then runs to left along the unstable manifold, as shown in Fig. 4(d), at last returns to the stable node, as shown in Fig. 4(b).

Two important characteristics can be found from Fig. 4(d). One is that the stable manifold of the saddle coincides exactly with the part of the threshold curve right to the stable node, and the other is that the trajectory running across the threshold curve results in an action potential. In addition, the unstable manifold (dashed green) drives the trajectory to move to right to form the action potential. Therefore, the dynamics of saddle, not only the stable manifold, but also the unstable manifold, play important roles to evoke a spike. The stable manifold of the saddle locating right to the stable node acts as the threshold curve to evoke an action potential from the stable node, which resembles that of the SNIC [16]. However, the threshold curve for the Hopf bifurcation has no direct relationship to the dynamics of the saddle. Therefore, the threshold curve for the BHom differs from that of Hopf bifurcation. Considering that the threshold curve is also different from that for the SNIC, the threshold curve for the BHom is novel.

The responses (upper panel) of the steady state to inhibitory pulse stimulations (lower panel) with strength $A$ = $-$0.5(blue), and $-$0.65(red), $-$0.66 (black), and $-$0.8 (green) are shown in Fig. 5(a), where the responses and inhibitory stimulations match with each other by the same color.

Different from the excitatory stimulation, the membrane potential decreases during the stimulation. The stronger the stimulation strength, the more negative the membrane potential at the termination time of the stimulation, as shown in Fig. 5(a). For weaker stimulations such as $A$ = $-$0.5 (blue) and $-$0.65 (red), after the stimulation, the membrane potential increases to a maximal value which is less than the membrane potential of the saddle $V_s$ (purple dashed line), then decreases and recovers to the resting potential. Not PIR spike but subthreshold membrane potential is induced by the weak inhibitory stimulation. For stronger stimulations such as $A$ = $-$0.66 (black) and $-$0.8 (green), the membrane potential can increase to be higher than $V_s$ to form an action potential, i.e. the PIR spike. After the spike, the membrane potential recovers to the steady state.

The trajectories corresponding to Fig. 5(a) ($A$ = $-$0.65 (red) and $-$0.66 (black)) are shown in Fig. 5(b). The red trajectory for subthreshold membrane potential ($A$ = $-$0.65) locates within both the yellow area and the black trajectory for the PIR spike ($A$ = $-$0.66) locates within the yellow and blanket areas. As shown in Fig. 5(c), the black trajectory during the stimulation runs across the left part of the threshold curve from up-right to down-left, whereas the red trajectory does not run across the threshold curve. The hollow square corresponds to the termination of the pulse stimulation. Therefore, in case of the inhibitory stimulation, it is through the down-left part of the threshold curve, which locates down-left to the stable node (red solid circle), that the suprathreshold trajectory transits from the subthreshold area to the suprathreshold area. The result shows that the down-left part of the threshold curve is responsible for the generation of the PIR spike. Such threshold curve for the PIR spike exhibits a negative slope (about $-$0.75) in ($V$, $w$) plane, as depicted in Fig. 5(c), and resembles that of the stable focus near Hopf bifurcation [16], where it is also the down-left part of the threshold sets that is responsible for the PIR spike.

After the stimulation, the red trajectory locates within the yellow area and very close to the threshold curve, approaches the left unstable manifold (dashed green line) of the saddle (half-filled circle) along the red arrow, then goes away from the saddle along the left unstable manifold, at last returns to the steady state (red solid circle). Different from the red trajectory, the black trajectory after the stimulation (square) locates within the blanket area and very close to the threshold curve first, as shown in Fig. 5(b) and (c), after approaching the saddle, exhibits dynamics resembling that of the excitatory stimulation to form action potential. The stable and unstable manifolds of the saddle still play roles in the formation of the PIR spike, as shown in Fig. 5(d).

3.1.3.   Threshold for coexisted stable node (right and close to BHom)

As shown in Fig. 1(a), when $u \in [u_{BHom}$, $u_{SN}]$, the FHN model exhibits the coexistence of stable node and stable limit cycle, which correspond to the steady state and spiking behavior, respectively. For example, Fig. 6 shows the coexisting behaviors at $u$ = $-$1.08 ($\in [u_{BHom}, u_{SN}]$), where Fig. 6(a) represents the changes of the membrane potential (upper) of the steady state and spiking behavior (lower) with respect to time, and Fig. 6(b) denotes the corresponding phase portraits. The blue line, red line, solid circle, half-filled circle, and hollow red circle are for nullclines $dV/dt$ = 0, $dw/dt$ = 0, the stable node, the saddle, and the unstable focus, respectively. The membrane potential of stable node and saddle is $V_n$ $\approx$ $-$0.96206499680548 and $V_s$ $\approx$ $-0.72853450846403$, respectively. The black curve in Fig. 6(b) represents the stable limit cycle corresponding to spiking behavior (lower) in Fig. 6(a).

As depicted in Fig. 6(b), the attraction domain of the stable node (red solid circle) is marked with yellow, and that of the stable limit cycle (black curve) is marked with blanket. That is to say, any initial value in yellow region induces FHN model (Eq. (1)) to generate subthreshold behavior, which converge to the stable node, and that in blanket region results in spiking behavior, which approaches the stable limit cycle. Therefore, the border between the yellow and white region has the role to separate subthreshold behavior and spiking behavior, which resembles that of mono-stable node shown in Fig. 3(b), where the border separates subthreshold behavior and suprathreshold (spike) behavior. The details of Fig. 6(b) around the saddle (half-filled circle) are enlarged in Fig. 6(c), where the solid (dashed) green line represents the stable (unstable) manifold of saddle (half-filled circle). There one striking characteristic is that the stable manifold (solid green line) of the saddle coincides with the border between the yellow and blanket regions, which resembles that shown in Fig. 4(d). Besides, the stable node shown in Fig. 6(b) is surrounded by the white region from left, below, and right side, which resembles to that for the mono-stable node shown in Fig. 3(b). These geometrical characteristics imply that the results for the coexisted node maybe resemble that of the mono-stable node.

For the stable node coexisting with spiking at $u$ = $-$1.08, excitatory suprathreshold stimulation such as $A$ = 0.23 (black) can induce the steady state changed to spiking, as shown in Fig. 7(a), while subthreshold stimulation such as $A$ = 0.22 (red) can not evoke the spiking behavior. The membrane potential of saddle $V_s$, denoted by purple dashed line, plays the role of voltage threshold curve for spiking. The corresponding phase trajectories and the attraction domain of the stable node (yellow) and spiking (white) are shown in Fig. 7(b), which exhibit dynamics similar to the mono-stable node shown in Fig. 4(b) to a large extent. The red trajectory for subthreshold stimulation ($A$ = 0.22) locates within the yellow area, and the black trajectory for suprathreshold stimulation ($A$ = 0.23) locates within both yellow and blanket areas. For $A$ = 0.23, the trajectory (black) during stimulation runs across the border between the attraction domain of the stable node (yellow)and the spiking (white), as shown in Fig. 7(c)), i.e. the stable manifold (solid green) of the saddle, as shown in Fig. 7(d). After the stimulation (hollow square), the trajectory (black) coincides with the stable limit cycle (thick magenta line) to form spiking behavior, as shown in Fig. 7(b). For $A$ = 0.22, the trajectory (red) during stimulation does not run across the stable manifold (bold green) of the saddle (half-filled circle), and after the stimulation recovers to stable node (red solid circle), as shown in Fig. 7(b), (c), and (d). Similar to the mono-stable node (Fig. 4), the stable manifold (solid green) of the saddle acts as threshold curve. Different from mono-stable node, spiking behavior instead of a spike is induced.

For the stable node coexisting with spiking at $u$ = $-$1.08, the inhibitory stimulation with $A$ = $-$0.52 (black) can induce the steady state (stable node) changed to spiking, with strength $A$ = $-$0.51 (red) cannot induced spiking but subthreshold behavior, as shown in Fig. 8(a). It can be clearly found that the membrane potential of saddle $V_s$, denoted by the purple dashed line, is the voltage threshold (at black arrow) for PIR spiking behavior, similar to Ref. [17]. The corresponding trajectories and the attraction domain of the stable node (yellow) and spiking (white) are shown Fig. 8(b), the enlargement of Fig. 8(b) around the phase point (square) at the termination time of the stimulation is illustrated in Fig. 8(c), and the enlargement of Fig. 8(b) near the saddle (half-filled circle) together with the stable (solid green) and unstable (dashed green) manifolds of the saddle are shown in Fig. 8(d).

As shown in Fig. 8(c), during the stimulation, the trajectory for suprathreshold stimulation ($A$ = $-$0.52) runs across the down-left part of the threshold curve, i.e., the border between the attraction domain of the stable node (yellow) and stable limit cycle (white), however, the trajectory for subthreshold stimulation ($A$ = $-$0.51) does not. Therefore, similarly to that of the mono-stable node shown in Fig. 5(d), the down-left part of the threshold curve, which is down-left to the stable node, is responsible for the PIR spiking, across which the inhibitory stimulation can induce spiking behavior and otherwise not. Besides, similar to mono-stable node shown in Fig. 5(d), the dynamics of the saddle are also important for the evoked spiking by inhibitory stimulation.

In addition, from the geometrical viewpoint, the threshold curve for spiking induced by excitatory stimulation (Fig. 7) and inhibitory stimulation (Fig. 8) can be attributed to that the stable node is surrounded all around by white area (the attraction domain of the stable limit cycle), which provides the chance that the stimulation drives the trajectory to transit from the yellow area (the attraction domain of the stable node) to the white area, i.e., to run across the border to form spiking. From viewpoint of the dynamics, it is the stimulation-induced transition from stable node to stable limit cycle.

3.2.   Threshold for stable focus near BHom

In this subsection, $c$ = $-$0.4, there still exists a BHom bifurcation with respect to parameter $u$. The resting state near BHom bifurcation now becomes stable focus. Two cases, the mono-stable focus (left and close to BHom) and the stable focus coexisting with spiking (right and close to Bhom) are considered. The results for the stable focus obtained in this subsection are quite similar to those for stable node depicted in subsection 3.1. Hence a lot of trivial details are omitted below for sake of simplicity.

3.2.1.   Bifurcations with respect to u and the BHom orbit

Shown in Fig. 9(a) is the bifurcation diagram with respect to $u$, which resembles Fig. 1(a) ($c$ = $-$0.55) except for the lower branch. The symbol $\rm {H_1}$ and $\rm{H_2}$ (black squares) at lower branch represent the supercritical Hopf bifurcation points, and the upper and lower green curves represent the extreme value of a subthreshold oscillation related to $\rm{H_1}$ and $\rm{H_2}$. Vmax (Vmin) represents the maximal (minimal) value of the stable limit cycle (red) bifurcated from $\rm{H_3}$, and stable limit cycle contacts with saddle to form a big homoclinic orbit (BHom) bifurcation, which occurs at $u_{BHom}$ $\approx$ $-$0.99447689769051 (black arrow) and is left to the $\rm{H_1}$ ($u_{H_{1}}$ $\approx$ $-$ 0.936). Fig. 9(b) represents the changes of firing period (blue) of the spike trains and the corresponding firing frequency (red) near BHom bifurcation (black arrow). Therein, it can be found that the value of firing period increases almost to infinity as $u$ decreased to $u_{BHom}$ (at black arrow). Meanwhile the firing frequency decreases to nearly zero. Therefore, resembling that shown in Fig. 1(b), such BHom bifurcation also exhibits type Ⅰ excitability.

The BHom orbit (black circle) at $u$ = $u_{BHom}$ for $c$ = $-$0.4 is illustrated in Fig. 10(a), which contains the stable focus (solid circle) and unstable focus (hollow circle) and goes exactly through the saddle (half-filled cycle). The red and blue lines are for the nullclines $dV/dt$ = 0 and $dw/dt$ = 0, respectively. The enlargement of Fig. 10(a) near saddle is depicted in Fig. 10(b), where the solid (dashed) green line represents the stable (unstable) manifold of saddle and the black arrow indicates the direction (anti-clockwise) for both the trajectory and manifolds of saddle, which further identifies that the homoclinic loop begins from and ends at the saddle indeed. The BHom orbit resembles that of $c$ = $-$0.55 except for containing a stable focus instead of stable node.

3.2.2.   Threshold for mono-stable focus (left and close to BHom)

When $u<u_{BHom}$, for example, when $u$ = $-$1.03 ($<u_{BHom}$), the system has mono-stability corresponding to stable focus. The membrane potential of the stable focus is $V_f$ $\approx$ $-$1.02368300429992. Suprathreshold excitatory stimulation such as $A$ = 0.58 (black) can induce an action potential followed by damping oscillations, while subthreshold excitatory stimulation such as $A$ = 0.57 (red) cannot induce an action potential but damping oscillations, as depicted in Fig. 11(a), which is different from the mono-stable node (Fig. 4), wherein the damping oscillation does not appear. Moreover, it is also the membrane potential of saddle $V_s$ ($\approx$ $-$0.436403782972816), denoted by purple dashed line, that acts as the voltage threshold for spike (at black arrow).

The phase trajectories of above two responses together with the threshold curve (i.e., the border between yellow and white region) are shown in Fig. 11(b). Here the yellow (white) region is for subthreshold (suprathreshold) area. The enlargement of dynamical behaviors between the stable focus (red solid circle) and saddle (half-filled circle) is depicted in Fig. 11(c), the details near the saddle (half-filled circle) are further enlarged in Fig. 11(d), where the hollow square is for the phase point at the termination time of stimulation. For $A$ = 0.58 (black), during the stimulation, the black trajectory runs across the threshold curve, i.e., the stable manifold of the saddle, to form a spike followed by damping oscillations, as shown in Fig. 11(c) and (d). However, for $A$ = 0.57 (red), the trajectory (red) cannot run across the threshold curve and thus forms the subthreshold damping oscillations. Similarly to the mono-stable node shown in Fig. 4, the stable manifold of the saddle acts the threshold curve for action potential evoked from the stable focus by excitatory stimulations. Different from the mono-stable node, the damping oscillations appear. Such difference is due to that the stable equilibrium in Fig. 4 is node, whereas here (Fig. 11) is focus.

In addition, although the result for the focus depicted in Fig. 11(d) in the present paper and Fig. 14 (right) in Ref. [16] seems similar with each other, they have difference in essence from the viewpoint of dynamics. The similarity is that, for focus near both BHom and Hopf bifurcation, the threshold curve for spike evoked by the excitatory stimulation locates right to the focus. However, the former in case of BHom bifurcation has close relationship with the saddle, i.e., coincides with the stable manifold of saddle, as shown in Fig. 11(d), whereas the latter in case of Hopf bifurcation has nothing to do with the saddle, as depicted in Fig. 14 (right) in Ref. [16]. Therefore, the threshold curve for the stable focus near the BHom bifurcation is a novel case different from that for the Hopf bifurcation.

The response of the resting state (stable focus) to inhibitory stimulations with strength $A$ = $-$0.59 (red) and $-$0.6 (black) are illustrated in Fig. 12(a), where the purple dashed line (i.e., the membrane potential $V_s$ of saddle) acts as the voltage threshold. As depicted at the black arrow, the evoked membrane potential for $A$ = $-$0.6 (black) increases to be higher than $V_s$ and then forms a single spike followed by the damping oscillations, whereas for $A$ = $-$0.59 (red) increases to nearly but fails to be higher than $V_s$ and then evolves to damping oscillations. The dynamical behaviors of the phase trajectories for $A$ = $-$0.59 (red) and $A$ = $-$0.6(black) are depicted in Fig. 12(b) together with subthreshold area (yellow) and suprathreshold area (white), where the squares denote the phase points at the end of stimulation. The two trajectories both start at and finally converge spirally to the stable focus (red solid circle) corresponding to the damping oscillation of the membrane potential shown in Fig. 12(a). The details at the squares are shown in Fig. 12(c), and around the saddle (half-filled circle) are shown in Fig. 12(d), where the solid (dashed) green line represents the stable (unstable) manifold of saddle. As depicted in Fig. 12(c), the trajectory for $A$ = $-$0.6 (black) runs across the threshold curve from up-right to down-left during the stimulation, whereas for $A$ = $-$0.59(red) does not. That is to say, the success or failure for the inhibitory stimulation to evoke a PIR spike depends on whether the trajectory during stimulation runs across the down-left part of threshold curve or not. Therefore, similar to that of mono-stable node in Fig. 5(d), it is the down-left part of the threshold curve, locating down-left to the stable focus, that acts as the threshold curve for PIR spike. Such result resembles that for focus near Hopf bifurcation [16], where it is also the down-left part of the threshold sets that are responsible for the PIR spike. Besides, as depicted in Fig. 12(d), one can also find that the manifolds of saddle also play important roles in the generation of the PIR spike in case for mono-stable focus.

3.2.3.   Threshold for coexisted stable focus (left and close to BHom)

As depicted in Fig. 9(a), when $u \in [u_{BHom}, u_{H_1}]$, the system exhibits the coexistence of stable focus and stable limit cycle, which correspond to the steady state and periodic spiking, respectively. For example, such coexistence at $u$ = $-$ 0.96 is depicted in Fig. 13. Shown in Fig. 13(a) is the resting potential (upper), which corresponds to the constant potential of stable focus, and period-1 spiking (lower), which corresponds to the stable limit cycle. Fig. 13(b) describes the corresponding phase portraits and its enlargement near saddle (half-filled circle) is depicted in Fig. 13(c). The intersections of the nullclines $dV/dt$ = 0 (red) and $dw/dt$ = 0 (blue) are stable focus (solid circle), saddle (half-filled circle), and unstable focus (hollow circle). The yellow (white) area in Fig. 13(b) is for the attraction domain of stable focus (limit cycle), and their border separates the subthreshold behavior and spiking behavior at parameter $u$ = $-$0.96, and the stable manifold (green solid) in Fig. 13(c) coincides exactly with the border near the saddle, which is quite similar to that of stable node coexisting with stable limit cycle shown in Fig. 6. The membrane potential of the saddle is $V_{s}$ $\approx$ $-$0.453557539050367 and that of the stable focus $V_{f}$ $\approx$ $-$0.953482960659727 at $u$ = $-$0.96.

In case of stable focus coexisting with the stable limit cycle at $u$ = $-$0.96, the excitatory stimulation can induce the resting state changed to spiking, as depicted in Fig. 14, which is different from that of the mono-stable focus shown in Fig. 11, where single spike is evoked. In Fig. 14(a), the red membrane potential ($A$ = 0.47) corresponds to the subthreshold potential and the black one ($A$ = 0.48) to the period-1 spiking, which implies that the constant membrane potential $V_s$, denoted by the purple dashed line, is the voltage threshold for the spiking behavior. Shown in Fig. 14(b) are the phase trajectories of above two responses together with the threshold curve, i.e., the border between the attraction domain of the stable focus (yellow) and stable limit cycle (white), where the thick magenta curve represents the stable limit cycle. The trajectories between the stable focus (solid circle) and saddle (half-filled circle) in Fig. 14(b) are further depicted in Fig. 14(c). The enlargement of Fig. 14(c) near saddle is shown in Fig. 14(d), together with the stable (solid green) and unstable (dashed green) manifolds of the saddle, where the stable manifold of saddle coincides with the threshold curve near saddle.

As can be found from Fig. 14(b)-(d), the black trajectory runs across the stable manifold and finally moves along the stable limit cycle (thick magenta line), which corresponds to spiking behavior, whereas the red trajectory does not run across the stable manifold and then forms subthreshold damping oscillation. Therefore, resembling the mono-stable focus shown in Fig. 11, the stable manifold of the saddle acts as the threshold curve for the spiking evoked by excitatory stimulation. Besides, the unstable manifold also plays important role in the generation of spiking. As shown in Fig. 14(d), the right part of the unstable manifold forces the black trajectory to move right and thus results in the increase of the membrane potential, which is a key factor to evoke spiking.

At $u$ = $-$0.96, the inhibitory stimulation with $A$ = $-$0.41 (black) induce PIR spiking, resembling Ref. [17], whereas with strength $A$ = $-$0.4 (red) cannot, as depicted in Fig. 15(a). The phase trajectories of above two responses are shown in Fig. 15(b) together with the threshold curve, i.e., the border between the attraction domain of stable focus (yellow) and stable limit cycle (white). Therein the hollow square denotes the phase point at the termination of stimulation and thick magenta line represents the stable limit cycle. The details around the hollow squares in Fig. 15(b) are shown in Fig. 15(c). There one observes that, during stimulation, the trajectory for $A$ = $-$0.41 (black) runs across the threshold curve and then spiking is formed, whereas for $A$ = $-$0.41 (red) does not and then no spiking is formed, which implies that the part of threshold curve down-left to the stable focus is responsible for the PIR spiking. The enlargement of Fig. 15(b) near the saddle are depicted in Fig. 15(d) combined with stable manifold (solid green) and unstable manifold (dashed green) of saddle. As can be found from Fig. 15(d), the manifolds of saddle also play important roles in the generation of spiking evoked by inhibitory stimulation in case for stable focus coexisting with the stable limit cycle.

3.3.   The two-parameter bifurcation diagram in ($u$, $c$) plane

In the present subsection, firstly, the bifurcations with respect to $u$ at different $c$ values, including BHom bifurcation, are acquired. Secondly, the two-parameter bifurcation diagram in ($u$, $c$) plane is acquired, which includes multiple codimension-1 bifurcation curves, such as BHom curve, and a codimension-2 bifurcation point related to BHom, i.e. the saddle-node Homoclinic orbit bifurcation point (SNHO). Finally, multiple dynamical behaviors separated by codimension-1 bifurcation curves in ($u$, $c$) plane, such as the parameter region for PIR spike and spiking, are acquired.

3.3.1.   Different bifurcations with respect to $u$ at different $c$ values

For the FHN model, except for the two cases of bifurcations with respect to $u$ (Fig. 1(a) for $c$ = $-$0.55 and Fig. 9(a) for $c$ = $-$0.4), two other cases of bifurcation with respect to $u$ at $c$ = $-$0.45 and $-$0.615 are depicted in Fig. 16(a) and (b), respectively. They have much similarity but small difference. Therefore, we shall focus on their difference, which mainly appears at the lower branch (LB) of $S$-shaped equilibrium curve. For $c$ = $-$0.45, the LB is formed by stable node (thick blue) and stable focus (thin blue), and a BHom bifurcation appears at the black arrow, as depicted in Fig. 16(a). At $c$ = $-$0.615, the whole LB corresponds to stable node, and the SNIC bifurcation (black arrow) occurs at the intersection of LB and MB, as shown in Fig. 16(b).

3.3.2.   The double-parameter bifurcations in ($u$, $c$) plane

To clearly show the 4 different cases of bifurcations with respect to $u$ at different $c$ values, the two-parameter bifurcation diagram in ($u$, $c$) plane is depicted in Fig. 17. There are five codimension-1 bifurcation curves, the curve for saddle-node bifurcation on an invariant cycle (SNIC, magenta line), curve for the saddle-node bifurcation (SN, blue), curve for big homoclinic bifurcation (BHom, black), curve for Hopf bifurcation on LB ($\rm{H_{LB}}$, dashed red line), and curve for Hopf bifurcation on UB ($\rm{H_{UB}}$, solid red). In addition, another important curve, shown by green line, represents the transition from stable node to stable focus on LB, and is short for NF. The intersection point between the curve SN and curve SNIC is a codimension-2 bifurcation point, saddle-node Homoclinic orbit bifurcation point (SNHO, cyan star), across which the curve BHom emerges. The SNHO has been reported in Refs. [16,18,34].

3.3.3.   Different dynamical behaviors in ($u$, $c$) plane

The ($u$, $c$) plane is divided into 7 regions by the above 6 curves, as shown by the different colors labeled with different roman numbers, the regions Ⅰ (orange), Ⅱ (gray), Ⅲ (green), Ⅳ (yellow), Ⅴ (white), Ⅵ (pink), and Ⅶ (cyan), as depicted in Fig. 17. Different regions correspond to different behaviors:

The region Ⅰ (orange) locates above the NF curve and left to BHom curve, which corresponds to mono-stable focus on LB.

The region Ⅱ (gray) lies below the NF curve, left to HBom curve and SNIC curve, which corresponds to mono-stable node.

The region Ⅲ (green) is upper to $\rm{H_{LB}}$ and corresponds to the coexistence of two stable limit cycles. One limit cycle exhibits small amplitude (Fig. 9(a)), which corresponds to subthreshold oscillation, and the other manifests large amplitude, which corresponds to spiking behavior.

The region Ⅳ (yellow) locates right to BHom curve, left to the right branch of the curve NF, and lower to curve the $\rm{H_{LB}}$, which corresponds to the coexistence of stable focus on LB and stable limit cycle (spiking behavior).

The region Ⅴ (white) is surrounded by the BHom curve, the NF curve, and SN curve, which corresponds to the coexistence of stable node and stable limit cycle.

The region Ⅵ (pink) locates right to the SN and SNIC, and left to $\rm{H_{UB}}$ curve, which corresponds to the stable limit cycle (spiking behavior).

The region Ⅶ (cyan) is right to the $\rm{H_{UB}}$ curve and corresponds to depolarization block.

In addition, the ($u$, $c$) plane is divided into four cases of bifurcations with respect to $u$ by three horizontal dashed lines, which go across codimension-2 bifurcation point (SNHO, cyan star) with $c$ $\approx$ $-$0.59, the minimal value ($c$ $\approx$ $-$0.548) of NF curve, and the minimal value ($c$ $\approx$ $-$0.419) of the $\rm{H_{LB}}$, respectively. From top to bottom, the four cases are labeled with case-1, case-2, case-3, and case-4, respectively. The bifurcations with respect to $u$ for $c$ = $-$0.55 in Fig. 1(a), and for $c$ = $-$0.4 in Fig. 9(a), and for $c$ = $-$0.45 in Fig. 16(a), and for $c$ = $-$0.615 in Fig. 16(b) correspond to case-1 to 4, respectively.

3.3.4.   The parameter region for the PIR spike in ($u$, $c$) plane

In the present paper, the mono-stable node and stable focus (close and left to the BHom curve) and the stable node or focus coexisting with stable limit cycle (close and right to the BHom curve) are investigated. The action potential, including PIR spike or spiking, can be evoked from these behaviors. The mono-stable equilibrium in the region marked by black star ($*$) can elicit PIR spike. The stable equilibrium in the region marked by red star ($*$) (i.e., Ⅳ and Ⅴ) can generate PIR spiking. For the stable equilibrium point coexisting with stable limit cycle and locating right to BHom curve, it is easy to understand the generation of the evoked spiking, which is the stimulation-induced transition from stable equilibrium to stable limit cycle. According to our results, it is natural to conclude that, in case of mono-stable equilibrium points, the action potential can be evoked, not just from the stable equilibrium near BHom, but also from all of those left to the BHom curve (black), as depicted in Fig. 17. The farther the stable equilibrium point left to the BHom, the stronger the critical strength of the stimulation to evoke action potential (verified already but not shown). However, for the stable node on the left of SNIC curve (magenta), as depicted in Fig. 17, it is thought to be difficult to evoke action potential by inhibitory stimulation at present.

4.   Conclusion and discussion

The nonlinear concept of threshold and post-inhibitory rebound spike are the fundamental conceptions in both nonlinear dynamics and neurophysiology, which are helpful to identify the physiological functions and modulation measures to firing pattern related to threshold or PIR spike [12,16,18,35]. In the present paper, we obtain the threshold curve for spike or spiking evoked from stable node and focus near bifurcation of big homoclinic (BHom) orbit in a modified FHN model without $I_h$ current, which confirms the conjecture in Ref [16]. The significances of the present paper are listed in the following aspects.

Firstly, an example of the threshold curve near the BHom bifurcation with type Ⅰ excitability is presented, as expected in Fig. 15 in Ref. [16]. For case of mono-stable node or focus close and left to the BHom bifurcation, the threshold curve is the border between the initial values in phase plane successful and unsuccessful to induce an action potential. For case of coexistence close and right to the BHom bifurcation, the threshold curve is the border between the attraction domains of the stable equilibrium and the stable limit cycle. For both the stable node and focus, the threshold curve is around the steady state from down-left, below, and right to the stable steady state, with structure resembling that of the focus near Hopf bifurcation to a certain extent. In future, the threshold curves for more kinds of bifurcations [16,18] should be acquired.

Secondly, the threshold curve for the BHom bifurcation is a novel case different from those of the well-known SNIC and Hopf bifurcation[16]. In the present paper, the part of the threshold curve right to the stable equilibrium coincides with the stable manifold of the saddle, which acts as the well-known threshold for spike or spiking evoked by excitatory stimulation. Such a result resembles that of the SNIC but is different from that of the Hopf bifurcations, wherein the threshold sets right to the stable focus is unrelated to saddle. For the BHom bifurcation, the part of the threshold curve down-left to the stable equilibrium, which has a negative slope, acts as the threshold for the spike or spiking induced by inhibitory stimulation, i.e. the PIR spike or spiking, which resembles that of the Hopf bifurcations. However, it is different from that of stable node near SNIC, whose threshold curve fails to extend to the region down-left to the stable node [16]. Therefore, the threshold curve for the BHom bifurcation presents a novel case of threshold curve. Although the threshold curves for many theoretical models are acquired [28], these curves have not been built relationship to bifurcations or PIR spike. In future, the characteristics of these threshold curves should be acquired and build relationships to bifurcations.

Thirdly, in the present paper, both parts of threshold curve down-left to and right to the steady state play important roles for the generation of PIR spike or spiking. The trajectory runs across the part of the threshold curve down-left to the steady state during the inhibitory stimulations, and then evolves to the neighborhood of saddle in the spike region. Near the saddle, the stable manifold of saddle, coinciding with the threshold, plays the role of separating the subthreshold behavior and subthreshold behavior. Furthermore, the unstable manifold plays the role to increase the membrane potential of the suprathreshold behavior by forcing the trajectory to move to right, which is a key step to form action potential or spiking. The detailed dynamical process and role for the threshold curve for the PIR spike are acquired.

Fourthly, the condition for PIR spike is extended from mainly the Hopf bifurcation with type Ⅱ excitability or $I_h$ current to BHom bifurcation unrelated to $I_h$ current. In most of the previous investigations, PIR spike in a neuron is related to either Hopf bifurcation or $I_h$ current [2,16]. Such a viewpoint is strengthened due to that $I_h$ current induces SN bifurcation changed to Hopf bifurcation point [35]. Furthermore, the condition for PIR spike is extended to be evoked from stable node near SNIC in system with $I_h$ current [12] and SN bifurcation near a Bogdanov-Takens bifurcation in system without $I_h$ current [18]. In the present paper, the condition for PIR spike is extended to be evoked from either node or focus near BHom bifurcation in a 2-dimensional neuron model without $I_h$ current. These results extend the condition for PIR spike to type Ⅰ excitability. In addition, the PIR spike or spiking are paradoxical phenomenon induced by inhibitory stimulations. Considering that the novel condition for PIR spike or spiking and multiple paradoxical phenomena induced by inhibitory effect [30] are helpful for identification of novel functions of inhibitory modulations in the nervous system [27], the relationship between paradoxical phenomenon such as the PIR spike and more kinds of bifurcations should be identified in future.

Last, the dynamical behavior near or related to the BHom bifurcation are acquired by double-parameter bifurcation analysis. In the double-parameter plane, multiple codimension-1 bifurcation curves, including the BHom curve, are acquired. The mono-stable steady state or steady state coexisting with spiking behavior are acquired, which is related to BHom bifurcation curve and PIR spikes. Furthermore, a codimension-2 bifurcation point related to the BHom bifurcation, the saddle-node Homoclinic orbit bifurcation point (SNHO) which is the intersection point between the codimension-1 bifurcation curves of SN, SNIC as well as BHom, is acquired. Therefore, the results of the present paper are for BHom bifurcation near the SNHO bifurcation. The SNHO bifurcation has been investigated in Refs. [4,16,18,20]. In addition, it is well-known that the common homoclinic orbit is related to another codimension-2 bifurcation, BT bifurcation [16,18,25,34]. In future, the relationship between the homoclinic orbit and SNHO or BT bifurcation should be studied.

The results of the present paper provide comprehensive viewpoint to the threshold for BHom, which extends the concept of threshold and the condition for PIR spike to a large extent, which are very important for both nonlinear science and neuroscience. In the present paper, the voltage threshold for a single pulse stimulation to evoke spike is investigated. In future, more stimulation parameters or stimulus patterns or current threshold such as stimulation directions, stimulation timings, stimulation amplitudes, and stimulation durations to evoke spike or spike patterns should be considered. In addition, the modulations of noise and time delay on the BHom orbit and the dynamics of network composed of neurons with BHom bifurcation should be studied.