Post inhibitory rebound spike related to nearly vertical nullcline for small homoclinic and saddle-node bifurcations

1.   Introduction

Nonlinear dynamics plays important roles in the identification of dynamical behaviors and functions of the nervous system ^[1,2]. Neuron is the basic unit of the nervous system and plays important roles in receiving, integrating, and transmitting information to achieve the corresponding biological functions via the action potential, or called spike for convenience ^[3,4]. In general, excitatory or positive stimulations enhance the membrane potential and thus are easy to evoke action potential, while the inhibitory or negative stimulations reduce the membrane potential and thereby are difficult to induce action potential ^[2,5]. However, it has been widely observed in the biological experiment that inhibitory stimulations induce some neurons with special currents to elicit action potential, i.e., the post-inhibitory rebound (PIR) spike ^[6,7]. In nonlinear dynamics, the PIR spike presents a typical example for the counterintuitive phenomenon in the nervous system, which has been build relationships to bifurcations or threshold curve for an action potential. In neuroscience, the PIR phenomenon is an intrinsic property of many neurons such as the CA1 pyramidal neurons ^[8] and deep cerebellar nuclei (DCN) neurons ^[9], and universally regarded as an important role in generating and maintaining biological rhythms such as feeding and swallowing behaviors ^[10]. Therefore, the study of PIR spike, including its biological functions, ionic current mechanisms, and dynamical mechanisms, has much practical importance in neuroscience and nonlinear dynamics.

More than a century ago, the excitable property of an inhibitory stimulation was observed in the response of a skeletal muscle to a reflex inhibition ^[11], which was identified by Sherrington (1913) who termed it "post-inhibitory excitation". Such a property was named as "post-inhibitory rebound" (PIR) by Granit (1956) ^[12], and also called rebound facilitation or anodal break response. The PIR phenomenon has been found to play important roles. For example, PIR phenomenon is related to locomotor rhythms in the spinal cord of rat ^[13] and snail feeding ^[14]. Besides, PIR phenomenon is also suggested to be an intrinsic "short-term memory" mechanism in the lateral pyloric neuron of crab ^[6] and to be related to sound rhythm coding in the Superior Paraolivary Nucleus in the auditory brainstem ^[7]. Moreover, it is reported that PIR phenomenon underlies a basis for theta rhythms ^[15] and depression ^[16]. Many underlying ionic mechanisms can account for the PIR phenomenon. For example, the hyperpolarization-activated cation current (I$_h$) is found to play important roles in generating PIR phenomena in various neurons such as the cerebral giant cells (CGCs) of gastropod mollusks ^[14] and individual brainstem neurons of mice ^[17]. There are other types of ionic currents, which underlie or contribute to the PIR phenomenon, such as the low-threshold T-type Ca$^{2+}$ current and A-type K$^{+}$ current ^[3]. In many cases, multiple mixed ionic currents are involved in the generation of PIR phenomenon via a complex process. In the leech HN cells, for example, hyperpolarization deinactivates low-threshold T-type Ca$^{2+}$ current and activates I$_h$, which produces a slowly depolarizing sag potential (one important characteristic of PIR phenomenon)^[18].

Except for the functions and ionic currents, the dynamical mechanisms for the PIR spike have also attracted much attention to present theoretical recognition to the PIR spikes ^[2,5,19]. In theory, the relationships between the PIR spike and codimension-1 bifurcations such as Hopf bifurcation (both subcritical and supercritical), saddle-node bifurcation on an invariant cycle (SNIC), and Big Homoclinic orbit (BHom) bifurcation have been acquired, and no obvious relationships to the Small Homoclinic orbit(SHom) bifurcation or saddle-node (SN) bifurcation have been acquired. In general, it has been widely accepted that the PIR spike is evoked from a stable focus near the Hopf bifurcation rather than a stable node ^[2,5]. For example, no PIR spike is evoked from a stable node near the SNIC bifurcation, based on the results of 2-dimensional Morris-Lecar model ^[5]. Recently, $I_h$ current induces SNIC bifurcation changed to Hopf bifurcation, which further implies that the PIR spike induced by $I_h$ current is related to the Hopf bifurcation ^[20]. However, more recently, in a three-dimensional neuronal model with $I_h$ current, the PIR spike appears near the SNIC bifurcation because the shape of threshold surface changes for the stronger $I_h$ current ^[21]. Moreover, PIR spike is conjectured to be associated with the BHom bifurcation in the well-known reference (Izhikevich, E.M., [2000]) ^[5], which then is numerically confirmed in the recent study ^[22]. The BHom orbit begins from and ends at a same saddle and is big enough to contain a steady state (stable node or focus). Another Homoclinic bifurcation involves a homoclinic orbit with small amplitude, which is called Small Homoclinic (SHom) orbit bifurcation for convenience. Different from BHom, the SHom orbit is relatively small and does not contain a stable steady state. Influenced by the well-known result that PIR spike can not be evoked from a node near the SNIC bifurcation in the ML model ^[5], it seems to be difficult to evoke a PIR spike from the stable node near the SHom or SN bifurcation. These results are helpful to develop comprehensive understandings of the PIR spike. However, up to now, no direct results of PIR spike have been provided for the SHom and SN bifurcations.

The generation of PIR spike has been well interpreted by the nonlinear concept, namely threshold, which is also an important characteristic of the action potential. For example, the well-known PIR spike, evoked from a stable focus near the Hopf bifurcation ^[5], can be attributed to the shape of the threshold curve. In phase plane with membrane potential as the horizontal axis, the threshold curve for an action potential surrounds the the stable focus (steady state) from left, below, and right sides (for more details, please see Figure 14 (right) in Ref ^[5]). Such shape of the threshold curve facilitates the generation of a PIR spike. An inhibitory stimulation (strong or prolong enough if necessary) can induce the trajectory initiating from the steady state to move to left and then run across the left part of the threshold curve, resulting in a PIR spike. Besides, in case of a stable node near SNIC bifurcation in the Morris-Lecar model ^[5], the cause for no PIR spike is that the threshold curve locates right to the stable node and has no a part locating down-left to the stable node, which is the cause for no PIR spike. The inhibitory stimulation-induced trajectory moves to left and thus has no chance to run across the threshold curve, resulting in no PIR spike (for more details, please see Figure 14 (left) in the study ^[5]). In fact, although not paid much attention, the nullcline of the recovery variable is nearly horizontal for no PIR spike whereas is non-horizontal for the PIR spike. The PIR spike for the BHom can be explained from the geometrical characteristic of the BHom orbit, i.e., containing the stable node. Such a characteristic provides the chance that the inhibitory stimulation-induced trajectory from the stable node runs across the down-left part of the BHom orbit, and thus is favorable to form a PIR spike. For the stable steady state prior to the BHom bifurcation, the threshold curve exhibits a shape similar to that of the stable focus prior to the Hopf bifurcation in Ref ^[5], which means that a PIR spike can be evoked. When talking about the SHom bifurcation, the geometrical characteristic of the SHom orbit (i.e., not containing a stable node) is unfavorable to form PIR spike. Similar result may be extended to the SN bifurcation due to that the stable node is also out of the coexisted stable limit cycle. However, the direct simulation results have not been presented for the threshold curve of stable node near the SHom as well as the SN bifurcation. Therefore, it remains unclear that the PIR spike can be evoked near the SHom bifurcation and SN bifurcation or not, which is an important question to be answered to present comprehensive recognitions to the PIR spike. If the answer is yes, then, can the generation of the PIR spike be interpreted with the threshold curve or other dynamical characteristics such as the vector fields and the shapes of nullcline (nearly horizontal or non-horizontal)? Furthermore, up to now, the generation of a PIR spike (or not) is acquired by numerical simulation and interpreted by analyzing the dynamical behaviors of the trajectory in phase space or phase plane ^[2,5,21,22], which is lack of analytical validations or explanations. Therefore, to present deep understandings for the dynamical mechanism of the PIR spike is very important for the neurodynamics.

In the present paper, we study the PIR spike and dynamical mechanism in the FitzHugh-Nagumo (FHN) model composed of two variables (One is the membrane potential $V$ and the other is the recovery variable $w$), due to that the nullcline of variable $w$ is non-horizontal. Firstly, for the stable node prior to the SHom bifurcation, PIR spike is simulated by a suitable inhibitory stimulation. Secondly, the PIR spike can be well explained by the dynamics of the threshold curve and vector fields. Furthermore, the PIR spike can be attributed to the nearly vertical characteristic of the lower part of $w$-nullcline in phase plane $(V, w)$. Our study suggests that such a characteristic of $w$-nullcline further induces dynamics of the vector fields favorable to form the PIR spike. Thirdly, based on the geometrical characteristic of $w$-nullcline, the membrane potential $V$ after the inhibitory stimulation is analytically estimated to be higher than the threshold of an action potential, which then implies the generation of a PIR spike. Besides, the parameter ranges for the PIR spike are also acquired. Last, the results of the PIR spike can be extended to the node near the SN bifurcation, which is close to the SHom bifurcation. Based on these results, the problems related to the PIR spike to be answered in future are discussed.

The rest of the present paper is organized as follows. Section 2 introduces the model, the concept of threshold curve, the inhibitory stimulation-induced response, and the method for integration as well as bifurcation. The results are depicted in Section 3 and includes the detailed dynamics of the PIR spike as well as the corresponding dynamical mechanisms. Section 4 displays the conclusion and discussion.

2.   Model and method

The FitzHugh-Nagumo (FHN) model was created by reducing the Hodgkin-Huxley model to two dimensions ^[23] and modified in Refs ^[2,5] as follows

where $V$ is for the membrane potential, $w$ for recovery variable, and $u$ is treated as the bifurcation parameter. Such an FHN model is dimensionless. The parameter $\epsilon$ = 1 and the function $s(w) = \frac{b}{1+e^{(c-w)/d}}$ is an $S$-shaped function, where $b$ = 2, $d$ = 0.05 and $c = -0.55$. The external stimulation $I_{stim}$ reads as

where $A$ ($< 0$) denotes the pulse strength. Such a pulse stimulation, stating at $t = 10$ and terminating at $t = 11$, is applied to the steady state to study the PIR spike. The pulse duration $1$ is chosen as a representative in the present paper, other values for the pulse duration or the starting time of the pulse do not change the results qualitatively.

In the present paper, if the membrane potential is greater than $1$, an action potential is thought to appear, otherwise no action potential but subthreshold potential. Then the phase plane $(V, w)$ is divided into two parts, which are marked by blank and yellow respectively in the present paper and respectively called spike part and subthreshold part for convenience. The spike part (blank) is for the collection of the initial value, initiating from which the trajectory forms an action potential, while the subthreshold part (yellow) for those which form no action potential but subthreshold potential. The border between these two parts is the so-called threshold curve, which was reported in Refs. ^[5,21,22].

Fourth-order Runge-Kutta method is utilized to perform the integration with time step 0.001. The bifurcation analysis is applied to Eqs (2.1) without pulse stimulation (i.e., $I_{stim} \equiv 0$) and performed by software XPPAUT ^[25].

3.   Result

3.1.   The bifurcation structure

The bifurcations with respect to $u$ of FHN model are depicted in Figure 1(a). Figure 1(b) is the enlargement of Figure 1(a). The blue $S$-shaped curve denotes the equilibrium points. The lower branch (LB) of the $S$-shaped curve is for stable node (bold solid blue), and the middle branch (MB) for saddle (dotted blue). The intersection of the MB and LB form the saddle-node (SN, black solid circle) bifurcation at $u_{SN}\approx -1.10632$, which denotes that the saddle on MB and stable node on LB intersect and then disappear with increasing $u$. The upper branch (UB) of the $S$-shaped curve represents the focus. As $u$ decreases, the stable focus (solid blue) changes to unstable focus (dashed line) via a supercritical Hopf bifurcation labeled with H (black hollow circle). Meanwhile a stable limit cycle appears, whose maximal and minimal values are represented by the upper and lower red curves respectively. As $u$ decreases, the stable limit cycle contacts with a saddle (black half-filled circle) on the MB to form the Small Homoclinic orbit (SHom) bifurcation at $u_{SHom}\approx -1.21253$ and then disappears. In addition, the insert figure in Figure 1(a) describes the local bifurcation diagram at the red arrow, which includes the stable node (solid blue line), the saddle-node bifurcation (red solid circle), the Hopf bifurcation (red hollow circle), and the SHom bifurcation (red half-filled circle), which are not the focus of the present paper.

Figure 2(a) represents the global view of the phase portrait at the SHom bifurcation. The $V$-nullcline (red, a shape-like cubic curve) intersects with the $w$-nullcline (blue) to form three equilibrium points, stable node (red solid circle), saddle (half-filled circle), and unstable focus (red hollow circle), which corresponds to the lower, middle, and upper branch of equilibrium curve shown in Figure 1, respectively. The cyan trajectory goes through the saddle (half-filled circle) and thus forms a homoclinic loop (cyan cycle).

The details for Figure 2(a) near saddle are enlarged in Figure 2(b), and plotted together with the manifolds of saddle, where the solid and dashed gray lines are for stable and unstable manifolds of saddle, respectively. The homoclinic orbit starts exactly from the saddle (half-filled circle), then moves to right and along the right part of unstable manifold (dashed gray), and finally returns to the saddle from the upper part of stable manifold (solid gray).

3.2.   Dynamics of threshold curve for stable node near SHom

The phase portraits of the resting state near the SHom bifurcation point is shown in Figure 3(a), where $u = -1.22$ is chosen as representative, which is left and close to $u_{SHom}\approx -1.21253$. The $w$-nullcline (blue) intersects with $V$-nullcline (red) to form stable node (red solid circle), saddle (half-filled circle), and unstable focus (red hollow circle). Furthermore, these two nullclines divide the phase plane into six regions, labeled as I, II, III, IV, V, and VI, respectively, which are introduced as follows:

Region I and V: upper to $w$-nullcline and lower to $V$-nullcline, $dV/dt > 0$ and $dw/dt < 0$, marked by lines paralleling to back diagonal.

Region II: lower to both $w$-nullcline and $V$-nullcline, $dV/dt > 0$ and $dw/dt > 0$, marked by lines paralleling to back diagonal.

Region III: lower to $w$-nullcline and upper to $V$-nullcline, $dV/dt < 0$, $dw/dt > 0$, marked by lines paralleling to back diagonal.

Region IV and VI: upper to both $w$-nullcline and $V$-nullcline, $dV/dt < 0$ and $dw/dt < 0$, marked by lines paralleling to diagonal.

For the symbols $(+, -)$, $(+, +)$, $(-, -)$ and $(-, +)$ in the phase plane, the former in brackets represents the sign of $dV/dt$, and the latter represents the sign of $dw/dt$. The simplified anti-clockwise rotating vector fields, formed by these dialog and anti-dialog lines in the phase plane, are related to the generation of action potential and are introduced in section 3.4. The enlargement around region III further depicted in Figure 3(b).

Figure 4(a) shows the threshold curve corresponding to the resting state at $u = -1.22$. If the phase point $(V, w)$ in the yellow area is chosen as an initial value, the behavior is subthreshold potential first and then recovers to the stable equilibrium point. If the initial phase point locates in blank area, an action potential appears at first and then recovers to the stable equilibrium point. The border between the yellow and blank areas is the so-called threshold curve for a spike or an action potential. Here the symbols are the same as those used in Figure 3.

The dynamics of threshold curve near the saddle are enlarged in Figure 4(b). The threshold curve runs across the saddle and coincides with the stable manifold (solid gray curve) of the saddle, which shows that the stable manifold (solid gray curve) of the saddle plays a role of the threshold curve. In addition, the unstable manifold (dashed gray curve) of the saddle coincides with $w$-nullcline (blue) near the saddle. The threshold curve here, i.e., the stable manifold of saddle, which locates right to the stable node, can be used to explain the well-known fact that an excitatory (positive or depolarization) stimulation can induce an action potential from the resting state, such as the stable node (red circle) in Figure 4(a). A suitable excitatory stimulation applied to the resting state induces the membrane potential $V$ to increase (i.e., from left to right), to run across the threshold curve near the saddle, and then to form an action potential. For convenience, $V = Vs$ (the membrane potential of saddle) can be taken as the voltage threshold near the saddle, which is a sufficient condition for the generation of a spike. The spike induced by an excitatory stimulation is very popular, but not the focus of the present paper.

The threshold curve down-left to the saddle is detailed in Figure 4(c), which extends into the area down-left to the stable node, not only into the region II but also into the region I. Then, it follows a natural consequence that, if the initial value of the model is chosen in the area down-left to the stable node and within the blank area of region I or II, an action potential appears. Based on the geometrical characteristic of the threshold curve, one can speculate about the generation of a PIR spike or not. When applied an inhibitory stimulation to the resting neuron (i.e., the stable node), the membrane potential decreases, that is to say, the phase point of the trajectory moves to left from the node in the phase plane. Notice that a part of the threshold curve locates left to the stable node. If the inhibitory stimulation, for example being strong or prolong enough, induces the behavior changed from the stable node to the blank area, then the trajectory runs within the blank area to form a PIR spike, and at last recovers to the stable node. Such speculation is evidenced in the following subsection.

3.3.   Dynamics for the PIR spike

In this subsection, we explore the dynamical mechanism of PIR spike from the perspective of threshold curve. Such a mechanism is obtained by studying the dynamics of trajectory in phase plane.

3.3.1.   The PIR spike and voltage threshold

For the stable node at $u = -1.22$ (corresponding to the stable node in Figure 3(a)), a PIR spike (cyan, upper panel) can be evoked by an inhibitory pulse stimulation with $A = -1$ (cyan, lower panel), as depicted in Figure 5. However, PIR spike can not be evoked by the the inhibitory pulse stimulation with $A = -0.99$ (red dotted). The insert figure in Figure 5 represents the enlargement around the arrow. Here, the stimulation strength $A = -1$ is chosen as the representative larger than the threshold of inhibitory current stimulation to evoke PIR spike, and $A = -0.99$ is chosen as representative lower than the threshold. The threshold strength to evoke a PIR spike is about $A\approx$ -0.999556, which is between $-1$ and $-0.99$.

As mentioned above, the membrane potential of saddle can be used as the voltage threshold of an action potential for convenience, and it does play such a role, as shown by the purple dashed line in Figure 5 ($Vs\approx -0.748796275688766$). During the stimulation, the membrane potential decreases. The stronger the stimulation, the more negative the membrane potential at the termination time of the stimulation, as shown in insert figure of Figure 5. After the stimulation, the responses for $A = -1$ (cyan) and $A = -0.99$ (red dotted) exhibit different dynamical behaviors. For $A = -0.99$ (red dotted), the membrane potential increases to a maximal value which is less than $Vs$, then decreases and recovers to the resting potential. For $A = -1$ (cyan), the membrane potential firstly increases to $Vs$, then further increases to be much larger than $Vs$ to form an action potential. After the action potential, the membrane potential decreases and recovers to the resting potential.

3.3.2.   Trajectory dynamics of the PIR spike

The phase trajectory (red dotted) for $A = -0.99$ is plotted with the threshold curve (i.e., the border of the yellow and blank area), $V$-nullcline (red) and $w$-nullcline (blue), and equilibrium points (red solid circle for the stable node and the half-filled circle for the saddle), as shown in Figure 6(a). Similarly, the phase trajectory (cyan) for $A = -1$ is shown in Figure 6(b). The hollow squares denote the phase points at the termination times of the stimulations. Both trajectories start at and finally recover to the stable node (red solid circle) in anti-clockwise direction, as shown by the arrows.

Figure 6(c) is the enlargement of two trajectories for $A = -1$ and $A = -0.99$ around the termination time (hollow square) of the stimulation. During the stimulation, the trajectory (cyan) for the PIR spike ($A = -1$) runs across the part of threshold curve locating down-left to the stable node, as shown in Figure 6(c), which is consistent with the speculation in subsection 3.2. Due to that the phase point labeled as cyan square locates within the blank area, the trajectory followed by the square runs from region I, to II, to VI, to IV to form the action potential, as shown in Figure 6(c) and (b). However, the trajectory for the subthreshold behavior ($A = -0.99$) does not run across the threshold curve, and thus results in that the whole trajectory (red dotted) is within the yellow area and forms subthreshold behavior instead of PIR spike, as shown in Figure 6(c) and (a). Therefore, from the geometrical viewpoint, the dynamical mechanism for the PIR spike is that, during the stimulation, the trajectory of the electrical behavior runs across the threshold curve in region I. Such threshold curve for the PIR spike resembles that of the focus near the Hopf bifurcation ^[5], where it is also the down-left part of the threshold curve that is responsible for the PIR spike.

Figure 6(d) represents the enlargement of both trajectories around the saddle (half-filled circle) and the stable and unstable manifolds of saddle (solid and dashed gray lines). The trajectory for $A = -0.99$ (red dotted) locates left to the threshold curve (i.e., the stable manifold of the saddle) and is within the yellow area. The trajectory for $A = -1$ (cyan) locates right to the threshold curve, and is within the blank area. After approaching the neighborhood of the saddle, the trajectory turns right to get far away from the saddle along the right part of the unstable manifold (dashed gray), resulting in a large increase of the membrane potential, which is the key factor to form a PIR spike. The results show that the membrane potential increases fast to form the action potential as the membrane potential is higher than the voltage of saddle ($Vs$), which is consistent with Figure 5.

Figure 6(e) and (f) shows the details of the phase trajectory in Figure 6(a) and (b) before returning to the stable node, respectively. Both trajectories run across region IV to recover to the stable node.

3.4.   Vector fields underlie the PIR spike

In this subsection, the vector fields are used to explain the dynamical mechanism of the PIR spike. Compared with threshold curve, vector fields show more details of the trajectory dynamics and provide more perceptual intuition for us to understand PIR spike. Shown in Figure 7(a) is the global view of the vector fields (back arrows), which are plotted with the nullclines, equilibrium points, and threshold curve. The details in the regions down-left to the saddle are further enlarged in Figure 7(b). Around the boundary between regions I and II, the vertical characteristic of the $w$-nullcline induces the horizontal distribution of the vector fields. These characteristics are favorable for the generation of PIR spike. For example, choose an initial value in region I, and then speculate how it flows according to the vector fields. Such a trajectory has the chance to flow into region III, even into the area right of the saddle (i.e., the generation of an action potential), which depends on the location of the initial value in region I. One can speculate that the lower the location of initial point in region I, the easier the generation of PIR spike.

Shown in Figure 8 are the vector fields (indicated by black arrows) plotted together with the phase trajectory for $A = -0.3$ (pink), $A = -0.8$ (black) and $A = -1.2$ (cyan). The squares represent the termination phases of the pulse stimulations. Via the hollow circles, the trajectories run across $w$-nullcline. The stimulations play the role to change the initial value from the stable node (red solid circle) to the phase points labeled as the squares, which then initiate different dynamical behaviors. The dynamics of these trajectories following the square, i.e., after the stimulation, are completely determined by the vector fields. It can be clearly found that the vector fields down-left to the stable node underlie three cases of trajectory dynamics, which are labeled as case-A, -B, and -C for convenience and described in details as follows:

(1) Case-A depicted by the pink trajectory. An inhibitory pulse stimulation with weak strength $A = -0.3$ induces the behavior changed from the stable node to a phase point (the square point). Then, the phase trajectory after the stimulation runs from the square point in region I, across the $w$-nullcline, to region II, and directly to stable node, resulting in subthreshold membrane potential instead of PIR spike.

(2) Case-B depicted by the black trajectory. After running across the $w$-nullcline, the phase trajectory firstly runs into region II, then into region III, and at last converges to the stable node. The response is still subthreshold membrane potential instead of PIR spike. The difference between case-B and case-A is the recovery process of the trajectory to the stable node. The trajectory enters into region III for case-B, whereas always runs in region II for Case-A.

(3) Case-C depicted by the cyan trajectory. The square for $A = -1.8$, i.e., the phase point at the termination time of the stimulation, now is within the blank area instead of yellow area. Different from case-A and -B, the trajectory for $A = -1.8$ runs across the threshold curve (in region I) from up-right to down-left during the pulse stimulation, which is consistent with the speculation in subsection 3.4. After the stimulation, the trajectory runs across the blank area of region I, $w$-nullcline, region II, and then reaches the region right to the saddle, runs along the unstable manifold of the saddle to form an action potential (in regions VI and IV, which is similar to Figure 6(b)), at last returns to the stable node. Compared to case-A and -B, one striking difference here is that the trajectory for case-C runs into the spike part (blank area) and thus the PIR spike is formed.

In short, weak stimulation fails to evoke PIR spike, such as case-A and -B, and strong stimulation succeeds in evoking PIR spike, such as case-C, where the vector fields are the underlying reason. The stronger the inhibitory stimulation, the more down-left the location of the square. The geometrical characteristic of the vector fields implies that the more down-left the location of the square, the larger amplitude the rotation of the square-initiating trajectory, then the easier the generation of the PIR spike. As matter of fact, these three cases of trajectories can be characterized by the location of the phase point on $w$-nullcline, as depicted by the hollow circles in Figure 8. The stronger the stimulation, the lower the location of the hollow circle, then the more favorable to form the PIR spike the vector fields. Once the hollow circle of a trajectory is low enough, for example, within the blank area, then the PIR spike appears. These results suggest that the lower part of $w$-nullcline plays an important role in the generation of PIR spike, which is analytically verified in next Subsection. The nearly vertical characteristic of $w$-nullcline is used in the proof process.

3.5.   The lower part of $w$-nullcline induces PIR spike

Suggested by the vector fields and dynamics of trajectory, the lower part of $w$-nullcline plays important roles in the generation of PIR spike. Compared with threshold curve and vector fields, the dynamics of $w$-nullcline is the essential cause for PIR spike. To better illustrate such a viewpoint, the inhibitory stimulation with strength $A = -2$ is chosen as a representative and the evoked trajectory is shown in cyan in Figure 9(a). The plots and symbols are similar with those used in Figure 6.

3.5.1.   The basic idea to prove the existence of PIR spike

After removing the inhibitory stimulation, the trajectory of electrical behavior runs across a point $E_{*}$ ($V_{*}, w_{*}$) locating on $w$-nullcline at time $t_{*}$, as shown by the cyan hollow circle in Figure 9(a). The PIR spike can be evoked when the $w_{*}$ value of point $E_{*}$ (determined by the stimulation strength $A$) is negative enough to a certain extent. To prove such an assertion, a pink area is chosen, as shown in Figure 9(b), which denotes the partial enlargement of Figure 9(a). For the pink area, the left boundary is chosen as $w$-nullcline, which corresponds to $E_{*}$; the right boundary is denoted by $V = V_{RB}$ (vertical dashed line) and its value is preset to be larger than membrane potential of the saddle (i.e., $Vs$), which corresponds to the voltage threshold for an action potential; the top boundary (upper dashed line) is denoted by $w = w_{TB}$, where $w_{TB} < -2/3$ is needed in the proof process (introduced in section 3.5.2); and the bottom boundary can be a value negative enough if necessary.

Our basic idea starts from considering a trajectory that initiates from $E_{*}$ with $w_{*}$ value negative enough. Then we prove that such a trajectory gets away from the pink area via not the top boundary but the right boundary. That is to say, the membrane potential of the trajectory can be higher than the voltage threshold of the action potential (i.e., the membrane potential of saddle $Vs$), which thus implies a PIR spike. Furthermore, a threshold for $w_{*}$, below which the corresponding trajectory surely runs across the right boundary of pink area (i.e., the appearance of PIR spike), can be acquired with our method. Our idea is realized through integration to an estimation described as follows:

which follows from the quotient of the two equations of FHN model without stimulation (i.e., Eqs (2.1) with $I_{stim} \equiv 0$).

The outline of the proof process is summarized as follows:

At first, in the pink area, one upper bound of $dw/dV$ is acquired by certain estimation and is labeled as $g(V)$, that is, $dw/dV\leq g(V)$ or $dw\leq g(V)dV$, where $g(V)$ is an integrable function.

Then, applying integration on both sides of $dw\leq g(V)dV$ from point $E_{*}$ ($V_{*}, w_{*}$) at time $t_*$ to a phase point $(V, w)$ at a certain time $t$ ($t > t_*$), one has $\int_{w_*}^{w}dw\leq \int_{V_*}^{V} g(V)dV$, which follows the estimation for $V$ and $w$ values, i.e., $w-w_* \leq G(V)-G(V_*)$, where the function $G(V)$ is the primitive function of $g(V)$. Based on such an inequality, the threshold and range of $w_*$ value to ensure the appearance of PIR spike are acquired at a fixed $u$ value ($u = -1.22$ as representative).

At last, for the steady states at different $u$ values, including those prior to the SHom and prior to the SN bifurcation shown in Figure 1, the threshold and range of $w_*$ value for the PIR spike are acquired.

3.5.2.   The upper bound $dw/dV\leq g(V)$

The integrable function $g(V)$ is acquired by the estimation to the upper bound of $dw/dV$ in Eqs (3.1), which is estimated by increasing its numerator and meanwhile decreasing its denominator in the pink area.

First, the range of $V$ value for the trajectory within the pink is estimated. The numerator of $\frac{dw}{dV}$, which is from $dw/dt = \epsilon (-u+V-s(w))$, is positive since the pink area locates in the region with $dw/dt > 0$ (region II). Therefore, $\epsilon (-u+V-s(w)) > 0$, where $\epsilon = 1$. Due to that $s(w) = \frac{b}{1+e^{(c-w)/d}}$ with $b$ = 2 is positive in the present paper, it can be acquired that $V > u+s(w) > u$. Note that the pink area locates left to the vertical dashed line $V = V_{RB}$, i.e., $V\leq V_{RB}$. Therefore, for the trajectory in the pink area, the membrane potential $V\in [u, V_{RB}]$.

Second, the lower bound of $dV/dt = V-V^3/3-w$ for the trajectory in the pink area is acquired. Let $f(V) = V-V^3/3$ and thus $df(V)/dV = 1-V^2$, which is zero for $V = 1$ and $V = -1$. $df(V)/dV$ is negative for $V\in [u, -1)$ and positive for $V\in (-1, V_{RB}]$, then, $f(-1) = -2/3$ denotes a local minimal value. In addition, note that $f(V)\geq -2/3 \Leftrightarrow V-V^3/3+2/3 \geq 0 \Leftrightarrow -\frac{1}{3}(V+1)^2(V-2)\geq 0 \Leftrightarrow V \geq 2$, where $\Leftrightarrow$ denotes the relation of equivalence. Then one has $\min\limits_{V\in [u, -V_{RB}]}f(V) = -\frac{2}{3}$ for any $V_{RB}\leq 2$. Due to that $w = w_{TB}$ is the top boundary of the pink area, then $dV/dt = f(V)-w\geq -2/3-w_{TB}$. To ensure $dV/dt > 0$, $w_{TB} < -2/3$ is needed. Therefore, one has a positive value of the lower bound of $dV/dt$, that is, $dV/dt\geq -2/3-w_{TB}$ for $V\in [u, V_{RB}]$ and $w \leq w_{TB}$, where $V_{RB}\leq 2$ and $w_{TB} < -2/3$ are needed.

Third, the positive value of the upper bound of $dw/dt = \epsilon (-u+V-s(w))$ for the trajectory in the pink area is acquired. Due to that $s(w) = \frac{b}{1+e^{(c-w)/d}}$ ($b = 2, c = -0.55$, and $d = 0.05$) exhibits monotonically increase with respect to $w$, then for any phase point on the cyan trajectory in pink area, one has $s(w)\geq s(w_*)$, which follows $dw/dt = \epsilon (-u+V-s(w))\leq \epsilon (-u+V-s(w_*))$. Besides, noticing that ($V_{*}, w_{*}$) is on $w$-nullcline ($dw/dt = 0$, i.e., $-u+V-s(w) = 0$), one has $s(w_{*}) = -u+V_{*}$. Thus, $dw/dt\leq \epsilon (-u+V-s(w_*)) = \epsilon (V-V_*)$.

It comes to show the upper bound of $dw/dV$, that is,

where $V, V_{*}\in [u, V_{RB}]$ and $V_{RB}\leq 2$, $w_{TB} < -2/3$. Therefore, the integrable function is acquired as

3.5.3.   The upper bound of the $w$ value of the trajectory in pink area

Integrate both sides of the last inequality from point $E_{*}$ ($V_{*}, w_{*}$) to a phase point $(V, w)$ within the pink area to form:

which then follows the upper bound of $w$, that is,

where the last "$\le$" is due to $V, V_{*}\in [u, V_{RB}]$. The result shows that $w$ of any phase point $(w, V)$ on the cyan trajectory in the pink area ($w\leq w_{TB}, V\leq V_{RB}$) is not higher than $-\frac{3 \epsilon (V_{RB}-u)^2}{2(2+3w_{TB})}+w_{*}$.

3.5.4.   The sufficient condition to ensure the generation of PIR spike

Note that the membrane potential of the trajectory is monotonically increasing in the pink area, which is due to $dV/dt > 0$ (proved already in section 3.5.2). Therefore, to make sure that the cyan trajectory runs across the right boundary of the pink area, it suffices to ensure that it has no chance to run across the top boundary $w_{TB}$, that is, to set

By equivalent transformation, such inequality can be rewritten into the following form

where $V_{RB}\leq 2$ and $w_{TB} < -2/3$. From the geometrical viewpoint, the solution $(V_{RB}, w_{TB})$ locates on or within an elliptic curve. Such an elliptic curve is centered at $V_{RB} = u$ and $w_{TB} = w_*/2-1/3$ with the length of horizontal semi-axis $-\frac{1}{2}\sqrt{\frac{2}{\epsilon}}(w_{*}+\frac{2}{3})$ and vertical semi-axis $-\frac{1}{2}(w_*+\frac{2}{3})$, which are both positive due to $w* \leq w_{TB} < -2/3$. A part of the elliptic curves for $u = -1.22$ and different $w_*$ values are depicted in Figure 10. Shown in orange, black, and green are the elliptic curves at $w_{*} = -1.3, -1.34,$ and $-2$, respectively. The vertical red dashed line corresponds to $Vs\approx-0.748796275688766$ and the hollow circle denotes the right vertex of each elliptic curve. For $w_{*} = -1.34$ (black) and $w_{*} = -2$ (green), the maximal value of $V_{RB}$ (hollow circle, denoted by $V_{RB\_max}$ below) is greater than $Vs$, which implies that the corresponding trajectory beginning from $E_*$ surely forms a PIR spike. More precisely to say, any $V_{RB}\in (Vs, V_{RB\_max}]$ can be chosen as the right boundary of the pink area in Figure 9(b) and used to prove the generation of PIR spike.

However, the result for $w_{*} = -1.3$ (orange in Figure 10) requires additional explanations. The maximal value of $V_{RB}$ occurs at orange hollow circle and is less than $Vs$, that is to say, $V_{RB}$ being greater than $Vs$ cannot be achieved. Therefore, our sufficient condition is not satisfied and thus no PIR spike can be proved with our developed method. It should be noticed that the point $E_*$ with $w_* = -1.3$ locates in the spike part (blank area), as depicted in Figure 9(b), and thus the trajectory initialing from $E_*$ is sure to form a spike with numerical simulation. The difference of the results between our developed method and numerical stimulation is explained in the following subsection 3.5.5.

The $V_{RB\_max}$ value obviously occurs at right vertex (hollow circle) of the elliptic curve and thus can be described as follows:

Therefore, for a fixed $w_{*}$, the inequality $V_{RB\_max} > Vs$ represents a sufficient condition to ensure the generation of PIR spike, and the interval $(Vs, V_{RB\_max}]$ provides the proper right boundary $V_{RB}$ in Figure 9(b) to be used in the proof process.

3.5.5.   The significance of the sufficient condition at $u=-1.22$

At $u = -1.22$, the sufficient condition now becomes

(1) The maximal $V_{RB}$ can be acquired to be negative correlation to the $w_{*}$, which shows that the more negative the value $w_{*}$, the larger the $V_{RB\_max}$. In fact, more strong external stimulation can induce more negative $w_{*}$. Therefore, the stronger the inhibitory stimulation, the easier the generation of the PIR spike. Such a result is consistent with the simulation results.

(2) The threshold of $w_{*}$ labeled as $w_{TH}$, below which the corresponding trajectory surely forms spike, is acquired. Due to $Vs \approx -0.748796275688766$, then $V_{RB\_max} > Vs$ is equivalent to $w_{*} < -1.333049364228327$. Therefore, based on our method, the threshold $w_{TH}$ is about $-1.333049364228327$. The result for $w_{*} = -1.34$ (black curve) depicted in Figure 10 corresponds to that near the threshold $w_{TH}$. It should be emphasized that such a threshold $w_{TH}$ obtained based on our developed method does not equal to the threshold acquired by the numerical simulation, which is about $w_{*}\approx-1.1$ (i.e., the lowest point of the yellow area depicted in Figure 9(b)). The threshold $w_{TH}$ of our method is lower than that of the numerical simulation, due to that the derivative of membrane potential is enlarged in our method.

(3) The more negative the value $w_{*}$, the longer the interval $(Vs, V_{RB\_max}]$, which implies that the wider range of the right boundary $V_{RB}$ can be used in the proof process.

3.5.6.   The PIR spike at different $u$ values

In fact, reviewing the proof process, for the steady state with $u$ value different from $u = -1.22$, the condition that

still ensures the generation of PIR spike if $Vs$ is chosen as the membrane potential of the saddle at $u$ (i.e., $V$ value of the middle branch of Figure 1(a)). The above inequality is equivalent to $w_{*} < \sqrt{2\epsilon}(u-Vs)-2/3$, which denotes the range of $w_{*}$ to evoke PIR spike at different $u$ values and is shown by the yellow area in Figure 11(a), where the black curve is for $w_{*} = \sqrt{2\epsilon}(u-Vs)-2/3$. In addition, the distribution of $V_{RB\_max}$ for the PIR spike at different values of $u$ and $w*$ are illustrated in Figure 11(b), where the color bar denotes the $V_{RB\_max}$ value, which is higher than the voltage threshold of spike, especially in the red area. Therefore, as long as $w_*$ value of the initial point $E_*$ is negative enough, the PIR spike is sure to appear at different values of $u$. Such a result holds true obviously for the steady state prior to not only the SHom bifurcation but also the SN bifurcation, as depicted by the two black arrows in Figure 11, which denote the positions of these two bifurcation points.

4.   Conclusions

The PIR spike is an important electrical behavior related to many biological functions such as rhythmic locomotor or information processing^{[1,2,3,4,10,13,14]}. The PIR spike is induced by the inhibitory stimulation, which is in contrast to the common view that a spike is always induced by an excitatory stimulation. Therefore, the PIR spike is an important counterintuitive nonlinear behavior to be studied in theory. In nonlinear dynamics, PIR spike has been related to the Hopf bifurcation or big homoclinic (BHom) orbit bifurcation^[2,5,22]. In the present paper, when the inhibitory stimulation is applied to the stable node near small homoclinic (SHom) orbit bifurcation and saddle-node (SN) bifurcation, the PIR spike is evoked. Such a PIR spike is further interpreted from the perspective of threshold curve, vector fields, and nullcline of recovery variable in the phase plane. The significances of the present paper are listed in the following aspects.

Firstly, the condition for PIR spike is extended to SHom and SN bifurcations. In the most of the previous investigations, PIR spike in a neuron is mainly related to either Hopf bifurcation or I$_h$ current ^{[2,5,14,17,21]}. Such a viewpoint is strengthened due to that I$_h$ current induces the SN bifurcation changed to the Hopf bifurcation ^[20]. Besides, for the resting state near the SNIC bifurcation, the stimulation fails to evoke PIR spike in a two-dimensional ML model ^[5], however succeeds in the three-dimensional modified I$_{Na, p}$+I$_K$ model with I$_h$ current ^[21]. The PIR spike is conjectured to appear for the BHom bifurcation ^[5] and confirmed in the FHN model ^[22]. Therefore, the result of the present paper extends the condition for the PIR spike to both SHom and SN bifurcations.

Secondly, the PIR spike is well interpreted by the nullcline of the recovery variable, threshold curve, and vector fields. In the existed literatures, threshold curve plays an important role to explain how the PIR spike is evoked and different shapes of threshold curve induce different PIR phenomena ^[5,22]. Such a viewpoint holds true in the present paper, which can be used to explain the PIR spike near both SHom and SN bifurcations. Except for the threshold curve, in the present paper, the dynamics of the vector fields and the nullcline of recovery variable are also used to account for the PIR spike. The key characteristic is that the lower part of the nullcline of recovery variable, which locates down-left to the stable node, is nearly vertical. Such a vertical characteristic induces that the dynamics of the vector fields are favorable for the PIR spike. The threshold curve and vector fields provide intuitive understanding of the PIR spike. The nullcline shows a new perspective to study the dynamics related to the PIR spike.

Last, different from the most of studies wherein the reason for the PIR spike is identified via numerical simulation, the role of nullcline played in the generation of PIR spike is approved analytically in the present paper. One sufficient condition for the PIR spike is acquired. Based on the parameter values, signs, and two functions in the FHN model, the equation to describe the derivative of membrane potential $V$ is estimated to acquire an integrable function. After integration to the integrable function, the range of $V$ after the inhibitory stimulation is acquired. The membrane potential $V$ can be larger than the voltage threshold for an action potential, then, PIR spike appears for the FHN model. Furthermore, the parameter ranges for the PIR spike are acquired. These results further provide insight into the dynamics of parameters which can be used to modulate the generation of the PIR spike.

In summary, the results of the present paper not only extend the conditions for the PIR spike to both SHom and SN bifurcations, but also reveal a new factor (i.e., the nullcline of the recovery variable) related to the PIR spike near the SHom and SN bifurcations, which are very important for both nonlinear science and neuroscience. Based on the results, in future, the problems related to the PIR spike to be answered are summarized in the following 4 aspects.

Firstly, the PIR spike can be evoked from high-dimensional models with the SHom or SN bifurcations (or not) should be further studied, which present comprehensive view to the emerging conditions for the PIR spikes or not.

Secondly, the influence of the nullcline on the generation of the PIR spike (or not) should be further studied, especially with different shapes of the nullclines for the SNIC bifurcation in the two-dimensional model. Up to now, it is widely accepted that PIR spike cannot be evoked from stable node near the SNIC bifurcation ^[5], where the nullcline of the recovery variable is almost horizontal. To the contrast, in the present paper, the vertical nullcline induces the generation of PIR spike. Therefore, in future, the shapes of nullcline should be regarded as an important factor to study the generation of PIR spike (or not) in 2-dimensional model with the SNIC bifurcation. The study of PIR spike should be accompanied by changing the shape of nullcline.

Thirdly, analytical proof process for the PIR spike (or not) near other bifurcation points such as the Hopf, SNIC, and BHom should be provided to present deep understandings. It is better to identify that the horizontal characteristic of the nullcline indeed induces no PIR spike for the SNIC bifurcation ^[5],

Last, the PIR spike in the network with the inhibitory coupling and time delay (or low decay rate) should be studied, due to that inhibitory coupling with slow decay rate (or time play) plays important roles in contributing to the synchronization of network such as Gamma rhythms ^[26].

Acknowledgements

This work is funded by the National Science Foundation of China (Grant Nos. 12072236 and 11872276).

Conflict of interest

The authors declare there is no conflict of interest.