Slow–Fast Dynamics Behaviors under the Comprehensive Effect of Rest Spike Bistability and Timescale Difference in a Filippov Slow–Fast Modified Chua’s Circuit Model

Abstract

Since the famous slow–fast dynamical system referred to as the Hodgkin–Huxley model was proposed to describe the threshold behaviors of neuronal axons, the study of various slow–fast dynamical behaviors and their generation mechanisms has remained a popular topic in modern nonlinear science. The primary purpose of this paper is to introduce a novel transition route induced by the comprehensive effect of special rest spike bistability and timescale difference rather than a common bifurcation via a modified Chua’s circuit model with an external low-frequency excitation. In this paper, we attempt to explain the dynamical mechanism behind this novel transition route through quantitative calculations and qualitative analyses of the nonsmooth dynamics on the discontinuity boundary. Our work shows that the whole system responses may tend to be various and complicated when this transition route is triggered, exhibiting rich slow–fast dynamics behaviors even with a very slight change in excitation frequency, which is described well by using Poincaré maps in numerical simulations.

1. Introduction

In the past decades, many research groups and scholars have studied slow–fast dynamical systems via experiments and computations based on extensive applications across many practical problems. For instance, neuronal activities and neural networks [1,2,3], calcium-ion oscillations in cells [4,5], catalytic oxidation reactions [6,7], and other applications [8,9,10].

Such systems frequently behave in complicated nonlinear phenomena characterized by various slow–fast dynamic behaviors such as relaxation oscillations, canard phenomena, bursting oscillations and mixed-mode oscillations. For more details, we suggest readers refer to the references [10,11] in addition to the references therein. In particular, bursting oscillations are generally more complicated and varied in oscillation patterns and generation mechanisms than the other three, and receive much more attention from scholars. In bursting oscillations, the waveform often presents periodic alternations between a cluster of spiking pattern (SP) and a period of quiescent state (QS). A general and powerful approach to deal with the mechanism of bursting phenomena in a slow–fast dynamical system is based on phase reduction and bifurcation theories [12,13]. Such methodology, i.e., the so-called slow–fast decomposition analysis, mainly consists of two steps: firstly, dividing the whole system into a fast subsystem and a slow subsystem via the timescale separation of state variables; secondly, utilizing the attractor structures and their bifurcations of the fast subsystem to explain the generation mechanism of potential bursting oscillations. Based on this method, bursting phenomena and the underlying generation mechanisms are well understood by employing dynamics of the fast subsystems. Relatively complete bursting patterns have been listed by Izhikevich [14] via slow–fast decomposition. More importantly, Izhikevich also provides a widely adopted classification method for bursting phenomena, by which one bursting oscillation can be named according to the two important bifurcations leading to transitions between QSs and the following SPs. For instance, “subcritical–Hopf/LPC” bursting means that the following dynamical generation mechanism can explain this bursting phenomenon: the trajectory will behave in SPs after it passes through a subcritical-Hopf bifurcation point and may enter into the following QS via the bifurcation LPC (fold bifurcation of limit cycles).

Since the dynamic features of the fast subsystems play the leading role in bursting oscillations, many researchers have paid increasing attention to the complicated structures existing in the fast subsystems, such as bistability and nonsmooth factors. Moreover, novel bursting patterns in addition to the related generation mechanisms have been explored in experiments and numerical computations in succession.

Bistability, referring to two coexisting attractors in a fixed parameter range, represents a common dynamics structure in nonlinear dynamical systems, and mounting evidence suggests that this geometric structure is a key factor (but not necessary) in the presence of bursting oscillations [2,14,15,16]. Two coexisting stable equilibria have been analytically investigated to be the principal reason of winner-take-all behavior in competition models [2,17], or the binocular rivalry in the so-called Levelt’s propositions [18]. Another bistability structure, namely rest spike bistability consisting of a stable limit cycle and a stable equilibrium separated by a saddle, has been studied as the core structure for slow spiking and bursting phenomena in the family of Morris–Lecar models [19]. The mixed bursting oscillation is a novel bursting form recently observed in a study by Ref. [20], where the trajectory can perform different SP patterns in a complete periodic motion. Duan et al. point out that such a special bursting pattern can be explained by bistability structures in the fast subsystem [21]. In continuous slow–fast dynamical systems, the slow–fast dynamic behaviors as well as the generation mechanisms are well understood based on slow–fast decomposition and the relatively well-developed conventional dynamics theories [22].

Discontinuity boundaries will exist in the fast subsystems where some nonsmooth events are involved, on which nonconventional dynamic behaviors may appear and participate in the dynamical evolution process of bursting oscillations. In nonlinear rotary drilling systems, friction induced by bit-rock interaction laws can cause the vector fields to become discontinuous, leading to the unique nonsmooth bursting patterns characterized by stick–slip motions [23,24,25]. Similar nonsmooth bursting patterns can also be observed when adopting a threshold control strategy in the Hindmarsh–Rose neuronal model [26]. On the other hand, owing to the unique advantages of nonsmoothness, almost all of the artificial neuron models are designed to be discontinuous [27].

In the slow–fast dynamical systems with discontinuous vector fields, complicated structures may be observed in the slow–fast dynamic behaviors when they interact with discontinuity boundaries, and the traditional dynamics and bifurcation theories may be useless to the understanding of those slow–fast dynamics behaviors. Fortunately, owing to the successful application of differential inclusion theory [28], nonconventional dynamics and bifurcation theories are developing rapidly (more details can be seen in [29] as well as the references therein). Riding on the nonconventional dynamics theories, scholars have proved that many nonconventional bifurcations may also affect the transition mechanisms of slow–fast dynamic behaviors, for instance, nonsmooth fold bifurcation and boundary homoclinic bifurcation can also cause the alternations between QSs and SPs [30,31,32]. Sliding bifurcations lead to unique sliding structures in SPs, although they do not result in transitions between QSs and SPs [33]. In [34], it was proved that boundary homoclinic bifurcation is a high codimension nonconventional bifurcation, which may result in different nonsmooth bursting patterns under different local structures of pseudo-equilibrium.

Almost all the above papers related to bursting phenomena indicate that bifurcations seem to be the only explanation for alternations between QSs and SPs. Recently, a different viewpoint was provided in the works [35,36,37] via proposal of the novel bursting mode referred to as pulse-shaped explosion, which can be attributed to the asymptotic behavior of vector field rather than bifurcations of attractors. Hastings et al. also pointed out that the saddle point and ghost attractor are the two essential elements leading to the non-bifurcation-induced transition in long-term ecological understanding [38]. This indicates that the related works regarding bursting oscillations also need to consider not only potential bifurcations but, more importantly, special flow structures in the fast subsystems.

In previous works related to bursting oscillations, the transition mechanism of bursting phenomena can be roughly classified into the following two types.

• The bifurcation-induced transition mechanism, which means that the bifurcation structures of attractors in the fast subsystems can explain the transitions in bursting oscillations very well.
• Non-bifurcation-induced transition mechanism, which means that the transitions in bursting oscillations can only be attributed to the unique manifold structures rather than to the bifurcation structures.

The bifurcation-induced transition mechanism in continuous (or discontinuous) slow–fast dynamical systems can be regarded as the replay of dynamics structures of bifurcations, which can be predicted well via bifurcation analyses of the fast subsystems, and the related achievements are plentiful and substantial. However, for the non-bifurcation-induced transition mechanism, there is not yet a useful research method, and the research also is relatively lag. Accordingly, further exploration of non-bifurcation-induced transition mechanism is important and meaningful to the comprehensive understanding of bursting phenomena in the slow–fast dynamics theory.

In the recent works by [30,39] involving Filippov slow–fast dynamical systems, unique nonsmooth rest spike bistability consisting of a stable nonsmooth cycle and the inner stable pseudo-equilibrium point is proved as the core structure of bursting patterns. Specifically, the alternation from the spike attractor (stable cycle) to the rest attractor (inner equilibrium) is not caused by bifurcations in the fast subsystems, thus indicating the occurrence of a non-bifurcation-induced transition. However, unfortunately, the underlying dynamics mechanism still needs to be clarified. Based on this, we will present further discussion of such non-bifurcation-induced transition structure in this paper by establishing a Filippov modified Chua’s circuit model with an external low–frequency excitation. Furthermore, we try to analyze the dynamics generation mechanism and the induced slow–fast dynamics behaviors, especially the influence on whole system responses, through quantitative calculations and qualitative analyses.

The rest of this paper is organized as follows. Section 2 presents a 3D Filippov slow–fast dynamical system by establishing a Filippov modified Chua’s circuit with an external low-frequency excitation. In Section 3, nonsmooth singularities are analyzed by employing nonconventional dynamics theory. Various slow–fast dynamics behaviors as well as the generation mechanisms are presented in Section 4 via numerical simulations. We draw some conclusions and discussions in Section 5.

2. Mathematical Model

The typical Chua’s circuit and modifications have been taken as examples in many studies of slow–fast dynamics via experiments or numerical simulations. To name but a few, canards and chaotic bursting were reproduced in memristor-based Chua’s circuit [40]; Marszalek and Trzaska analyzed mixed-mode oscillations and relaxations in slow–fast modified Chua’s circuit models with autonomous vector fields [41]. Particularly, those modified Chua’s circuits with low–frequency excitations are also employed to investigate various bursting oscillations and the generation mechanisms [24,25,26,42].

Here, we consider one modified Chua’s circuit by introducing a current source, as shown in Figure 1a. Such a circuit system consists of two capacitances $C_{1,2}$, one inductance L, one constant resistance R, one sinusoidal current source $i_S=I_{m}sin\left(\omega t\right)$ and one nonlinear resistance $R_N$. According to the Kirchhoff Law, the state equations can be expressed in the 3D ordinary differential equations

$$\left\{\begin{array}{c}{C}_1\frac{\mathrm{d}{u}_1}{\mathrm{d}\mathrm{t}}=\frac{u_2-u_1}{R}-g\left(u_1\right)+i_S,\hfill \\ C_2\frac{\mathrm{d}{u}_2}{\mathrm{d}\mathrm{t}}=i_L+\frac{u_1-u_2}{R},\hfill \\ L\frac{\mathrm{d}{i}_L}{\mathrm{d}\mathrm{t}}=-u_2,\hfill \end{array}\right.$$

(1)

where $R_N$ is designed as a nonlinear resistance possessing a piece-wise smooth voltage–current characteristic via introducing one comparator, as shown in Figure 1b. Then, the voltage–current characteristic can be written as

$$\begin{array}{cc}\hfill i_N=g\left(u_1\right)=& \hfill \frac{R_{12}(R_{15}+R_{16})+R_{13}{R}_{16}}{R_{12}{R}_{15}{R}_{18}}{u}_1-\frac{R_{10}{R}_{16}(R_7+R_6)}{100R_6{R}_8{R}_{11}{R}_{18}}{u}_1^3+\frac{(R_1+R_2)R_{16}}{R_1{R}_4{R}_{18}}{E}_{sat}\mathrm{sgn}\left(u_1\right).\hfill \end{array}$$

(2)

It can be found that $u_1=0$ defines a discontinuity boundary, labeled as $\mathrm{\Sigma }$, which means that the vector field of system (1) then is piece-wise smooth.

Using the following transformations $\tau =\frac{t}{R{C}_2}$, $x=\frac{u_1}{E_{sat}}$, $y=\frac{u_2}{E_{sat}}$, $z=\frac{R{i}_L}{E_{sat}}$, $\alpha =\frac{C_2}{C_1}$, $a=\frac{(R_{12}{R}_{15}+R_{12}{R}_{16}+R_{13}{R}_{16})R}{R_{12}{R}_{15}{R}_{18}}$, $b=-\frac{R_{10}{R}_{16}(R_7+R_6)R{E}_{sat}^2}{100R_6{R}_8{R}_{11}{R}_{18}}$, $\delta =-\alpha \frac{(R_1+R_2)R_{16}R}{R_1{R}_4{R}_{18}}$, $A=\frac{R{I}_m\alpha }{E_{sat}}$, $\Omega =\omega R{C}_2$ and $\eta =\frac{R^2{C}_2}{L}$, we can rewrite system (1) into the dimensionless form

$$\left\{\begin{array}{c}\dot{x}=\alpha \left(y-x-f\left(x\right)\right)+\delta \mathrm{sgn}\left(H\left(\mathbf{X}\right)\right)+Asin(\Omega \tau ),\hfill \\ \dot{y}=x-y+z,\hfill \\ \dot{z}=-\eta y,\hfill \end{array}\right.$$

(3)

where $·=\frac{d}{d\tau }$, $f\left(x\right)=ax+b{x}^3$, $\mathbf{X}={\left[xyz\right]}^{\mathbf{T}}$, $\alpha$, b, A, $\Omega$, $\eta$ are positive while a, $\delta$ are negative, $H\left(\mathbf{X}\right):=x=0$ represents the discontinuity boundary $\mathrm{\Sigma }$.

With the assumption that $\Omega$ is small enough, i.e., $0<\Omega \ll 1$, then there exists an order gap between the rates of traditional variables (x, y and z) and the whole excitation, admitting that the system (3) may perform slow–fast dynamics behaviors. However, there is no apparent fast variables and slow variables in (3), indicating that the classical slow–fast decomposition method cannot be directly applied to (3) to explain the potential slow–fast dynamics behaviors such as bursting oscillations since the division of the slow subsystem and the fast subsystem in (3) is the first requirement for the application of this method. In order to apply the classical slow–fast decomposition method to deal with the potential slow–fast dynamics behaviors in (3), we now introduce the following rise–dimension method to realize the timescale separation in (3).

Actually, $\left(Asin(\Omega \tau ),Acos(\Omega \tau )\right)=\left(w\left(\tau \right),v\left(\tau \right)\right)$ can be described by a simple planar–linear system ($\dot{w}=\Omega v;\dot{v}=-\Omega w$) with the initial condition $\left(w\left(0\right)=Asin\left(0\right),v\left(0\right)=Acos\left(0\right)\right)=(0,A)$. Therefrom, system (3) then is equivalent to a 3D Filippov fast subsystem (FS) controlled by a 2D smooth slow subsystem (SS), respectively, expressed as

$$\left\{\begin{array}{c}\dot{x}=\alpha \left(y-x-f\left(x\right)\right)+\delta \mathrm{sgn}\left(H\left(\mathbf{X}\right)\right)+w,\hfill \\ \dot{y}=x-y+z,\hfill \\ \dot{z}=-\eta y,\hfill \end{array}\right.\left(\mathrm{FS}\right)$$

(4)

and

$$\left\{\begin{array}{c}\dot{w}=\Omega v,\hfill \\ \dot{v}=-\Omega w.\hfill \end{array}\right.\left(\mathrm{SS}\right)$$

(5)

We then have transformed (3) into the standard form represented by a 5D autonomous Filippov slow–fast dynamical system that embodies three fast variables (x, y and z) and two slow variables (w and v). Although the rise–dimension method seems to increase the difficulty in the analysis of (3), it should be pointed out that only the slow variable w appears in the fast subsystem (4), indicating that the dynamics structures of the fast subsystem only can be affected by w. As a consequence, the classical slow–fast decomposition method can be directly employed to explain the generation mechanism of slow–fast dynamics behaviors in (3) just via regarding the slow variable w as a bifurcation parameter in (4). Meanwhile, the rise–dimension method also demonstrates that the generalized slow–fast analysis method proposed in [42] is feasible mathematically and logically.

When attractors of the fast subsystem come into contact with the discontinuity boundary $\mathrm{\Sigma }$, unique nonsmooth dynamics such as sliding bifurcations and boundary equilibrium bifurcations may appear and affect the responses of whole system (3). In the following section, we mainly focus on the analyses of nonsmooth dynamics on $\mathrm{\Sigma }$.

3. Nonsmooth Dynamics Analyses

Note that the discontinuity boundary $\mathrm{\Sigma }$ divides the whole phase space $R=\left\{\mathbf{X}\right|(x,y,z)\}$ into two subregions $R_{+}=\left\{\mathbf{X}\right|(x,y,z),(x>0)\}$ and $R_{-}=\left\{\mathbf{X}\right|(x,y,z),(x<0)\}$, where two 3D smooth subsystems, denoted by $S_{+}(x>0)$ and $S_{-}(x<0)$, can be observed in the two subregions, and the vector fields are, respectively, labeled as ${\mathbf{F}}_{+}(x>0)$ and ${\mathbf{F}}_{-}(x<0)$. Let ${\mathbf{X}}_{\mathrm{\Sigma }}(0,y_{\mathrm{\Sigma }},z_{\mathrm{\Sigma }};w)$ be a point on $\mathrm{\Sigma }$, by using the two Lie derivatives ${\mathcal{L}}_{{\mathbf{F}}_{+}}H=\langle \nabla H,{\mathbf{F}}_{+}\rangle$ and ${\mathcal{L}}_{{\mathbf{F}}_{-}}H=\langle \nabla H,{\mathbf{F}}_{-}\rangle$ ($\nabla H=\frac{\partial H\left({\mathbf{X}}_{\mathrm{\Sigma }}\right)}{\partial {\mathbf{X}}_{\mathrm{\Sigma }}}$, $\langle ·,·\rangle$ denotes the vector inner product), the discontinuity boundary $\mathrm{\Sigma }$ can be divided into different regions [43].

If ${\mathcal{L}}_{{\mathbf{F}}_{+}}H·{\mathcal{L}}_{{\mathbf{F}}_{-}}H>0$, i.e., ${\mathbf{F}}_{+}$ and ${\mathbf{F}}_{-}$ have the same direction on $\mathrm{\Sigma }$ at ${\mathbf{X}}_{\mathrm{\Sigma }}$, then we can define the two so-called sewing regions

$$\left\{\begin{array}{c}\hfill {\mathrm{\Sigma }}_{+}=\left\{{\mathbf{X}}_{\mathrm{\Sigma }}\right|\alpha y_{\mathrm{\Sigma }}+w-\delta >0,\alpha y_{\mathrm{\Sigma }}+w+\delta >0\},\hfill \\ \hfill {\mathrm{\Sigma }}_{-}=\left\{{\mathbf{X}}_{\mathrm{\Sigma }}\right|\alpha y_{\mathrm{\Sigma }}+w-\delta <0,\alpha y_{\mathrm{\Sigma }}+w+\delta <0\},\hfill \end{array}\right.$$

(6)

where the trajectory will cross through the discontinuity boundary without sliding motion when it interacts with $\mathrm{\Sigma }$ in ${\mathrm{\Sigma }}_{+}$ (or ${\mathrm{\Sigma }}_{-}$).

If ${\mathcal{L}}_{{\mathbf{F}}_{+}}H·{\mathcal{L}}_{{\mathbf{F}}_{-}}H<0$, ${\mathcal{L}}_{{\mathbf{F}}_{+}}H<0$ and ${\mathcal{L}}_{{\mathbf{F}}_{-}}H>0$, i.e., both ${\mathbf{F}}_{+}$ and ${\mathbf{F}}_{-}$ point to $\mathrm{\Sigma }$ at ${\mathbf{X}}_{\mathrm{\Sigma }}$, then we can define a so-called sliding region

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_S=\left\{{\mathbf{X}}_{\mathrm{\Sigma }}\right|\alpha y_{\mathrm{\Sigma }}+w-\delta >0,\alpha y_{\mathrm{\Sigma }}+w+\delta <0\},\end{array}$$

(7)

where the trajectory will enter into $\mathrm{\Sigma }$ and behave in the special sliding motion when it interacts with ${\mathrm{\Sigma }}_S$.

In particular, there are two so-called sliding boundaries between the two sewing regions (6) and the sliding region (7), which can be defined by ${\mathcal{L}}_{{\mathbf{F}}_{+}}H=0$ and ${\mathcal{L}}_{{\mathbf{F}}_{-}}H=0$, respectively, expressed by

$$\left\{\begin{array}{c}\hfill \partial {\mathrm{\Sigma }}_{S+}=\left\{{\mathbf{X}}_{\mathrm{\Sigma }}\right|\alpha y_{\mathrm{\Sigma }}+w+\delta =0\},\hfill \\ \hfill \partial {\mathrm{\Sigma }}_{S-}=\left\{{\mathbf{X}}_{\mathrm{\Sigma }}\right|\alpha y_{\mathrm{\Sigma }}+w-\delta =0\}.\hfill \end{array}\right.$$

(8)

Meanwhile, one 2D sliding vector field located in ${\mathrm{\Sigma }}_S$ can be derived as ${\mathbf{F}}_S=\frac{{\mathbf{F}}_{-}+{\mathbf{F}}_{+}}{2}+Q\frac{{\mathbf{F}}_{+}-{\mathbf{F}}_{-}}{2}$ by introducing a scalar $Q\in [-1,1]$ via the Utkin’s equivalent control method, corresponding to a linear-planar system

$$\left\{\begin{array}{c}{\dot{y}}_{\mathrm{\Sigma }}=-y_{\mathrm{\Sigma }}+z_{\mathrm{\Sigma }},\hfill \\ {\dot{z}}_{\mathrm{\Sigma }}=-\eta y_{\mathrm{\Sigma }}.\hfill \end{array}\right.$$

(9)

Besides the two 3D smooth subsystems $S_{+}$ and $S_{-}$, one may find that the Filippov fast subsystem (4) also has one 2D sliding subsystem (9) in ${\mathrm{\Sigma }}_S$. Modulated by the slow variable w, attractors such as equilibria and cycles of smooth subsystems may collide with a sliding boundary in (8), behaving in nonconventional bifurcations. In other words, nonconventional bifurcations in (4) generally occur on sliding boundaries. We now turn to discuss the equilibria and their stabilities as well as nonconventional bifurcations.

3.1. Equilibria and the Stabilities

Since the system (3) is $Z_2$ symmetrical, the equilibrium points of the fast subsystem (4) can be written as ${\mathbf{X}}_0=(x_0,0,-x_0;w)$. The equilibria in Filippov dynamical systems can be classified into the following four types: admissible equilibrium, pseudo equilibrium, boundary equilibrium and virtual equilibrium [29]. Furthermore, the first three types of equilibria really exist in the state space. In contrast, the last type of equilibrium does not exist in the state space and should be useless for understanding the slow–fast dynamics behaviors. Based on the fact that the Filippov fast subsystem (4) consists of two 3D smooth subsystems $S_{+}$, $S_{-}$ as well as one 2D sliding subsystem, we then have the following definition for the equilibria that really exist in the fast subsystem (4).

Definition 1.

The real equilibria of the fast subsystem (4) can be classified into the three types:

1. 

if ${\mathit{F}}_{+}\left({\mathit{X}}_0\right)=\mathit{0}$ and $x_0>0$( ${\mathit{F}}_{-}\left({\mathit{X}}_0\right)=\mathit{0}$ and $x_0<0$), ${\mathit{X}}_0$ is a so-called admissible equilibrium (AE), and $x_0$ is restricted by $w=x_0+f\left(x_0\right)-\delta$($w=x_0+f\left(x_0\right)+\delta$);

2. 

if $\left|w\right|<\left|\delta \right|$ and $x_0=0$, then ${\mathit{X}}_0$ is a so-called pseudo equilibrium (PE) located in sliding region ${\mathrm{\Sigma }}_S$, and which is equal to the equilibrium of sliding subsystem (9);

3. 

if ${\mathit{F}}_{+}\left({\mathit{X}}_0\right)=\mathit{0}$, $x_0=0$ and $w=-\delta$(${\mathit{F}}_{-}\left({\mathit{X}}_0\right)=\mathit{0}$, $x_0=0$ and $w=\delta$), ${\mathit{X}}_0$ is a so-called boundary equilibrium (BE) located on $\partial {\mathrm{\Sigma }}_{S+}$ ($\partial {\mathrm{\Sigma }}_{S-}$).

Note that AEs (labeled as $A{E}_{\pm }$) and PEs are the equilibrium points of the 3D smooth systems $S_{\pm }$ and the 2D smooth system (9), respectively, the corresponding characteristic equations of which can determine the stabilities according to the Routh–Hurwitz Criterion.

For an AE, the corresponding characteristic equation can be computed as the unified form

$$f_{AE}\left(\lambda \right):={\lambda }^3+(\kappa +2){\lambda }^2+(\kappa +\eta )\lambda +\eta (\kappa +1)=0,$$

(10)

where $\kappa =\alpha \left(a+3b{x}_0^2\right)$. Then we have the following result: if $(\kappa +2)>0,(\kappa +2)(\kappa +\eta )-\eta (\kappa +1)>0$ and $\eta (\kappa +1)>0$ are satisfied, the AE is stable.

Meanwhile, the stability of PE is characterized by the corresponding characteristic equation

$$f_{PE}\left(\lambda \right):={\lambda }^2+\lambda +\eta =0.$$

(11)

Then the two eigenvalues are easy to be derived as ${\lambda }_{1,2}=-0.5\pm 0.5\sqrt{1-4\eta }$, and PE is always stable when $\eta >0$. Furthermore, it is a stable node when $0<\eta \le 0.25$, while it is a stable focus when $\eta >0.25$.

3.2. Boundary Equilibrium Bifurcations

The two boundary equilibria, computed as ${\mathbf{X}}_{BE\pm }=(0,0,0;\mp \delta )$, will connect the AE and the PE when $A>-\delta$, at which codimension-1 boundary equilibrium bifurcations (BEBs) occurs. Taking the boundary equilibrium ${\mathbf{X}}_{BE+}$ as an example, codimension-1 BEBs need to meet the two following conditions [34]

$$\left\{\begin{array}{cc}\hfill \det \left({\mathbf{F}}_{+,\mathbf{X}}\left({\mathbf{X}}_{\mathrm{BE}+}\right)\right)=& \hfill -\eta \alpha (1+a)\ne 0,\hfill \\ \hfill \frac{\mathrm{d}H\left({\mathbf{X}}_{\mathrm{BE}+}\right)}{\mathrm{d}w}=& \hfill \frac{1}{\alpha (1+a)}\ne 0,\hfill \end{array}\right.$$

(12)

where the first one and the second one can be regarded as the nondegeneracy condition and the transversality condition, respectively.

In the arbitrary small neighborhood of ${\mathbf{X}}_{BE+}$, there exist an AE (labeled as ${\mathbf{X}}_{AE+}=(x_0,0,-x_0;w_1)$$(x_0>0)$) and a PE (labeled as ${\mathbf{X}}_{PE}=(0,0,0;w_2)$). Obviously, we then have

$$\left\{\begin{array}{c}{\mathbf{F}}_{+}\left({\mathbf{X}}_{AE+}\right)=\mathit{0},\hfill \\ H\left({\mathbf{X}}_{AE+}\right):=x_0>0,\hfill \end{array}\right.$$

(13)

at ${\mathbf{X}}_{AE+}$ and

$$\left\{\begin{array}{c}{\mathbf{F}}_{+}\left({\mathbf{X}}_{PE}\right)+q({\mathbf{F}}_{-}\left({\mathbf{X}}_{PE}\right)-{\mathbf{F}}_{+}\left({\mathbf{X}}_{PE}\right))=\mathit{0},\hfill \\ H\left({\mathbf{X}}_{PE}\right):=0,\hfill \end{array}\right.$$

(14)

at ${\mathbf{X}}_{PE}$ by introducing the transformation $q=\frac{1-Q}{2}$, respectively. Linearizing (13) and (14) at ${\mathbf{X}}_{BE+}$, we can find that

$$q=-\frac{x_0\alpha \left(1+a\right)}{2\delta }.$$

(15)

Noting that $x_0>0$, we have the following definition for codimension-1 BEBs at ${\mathbf{X}}_{BE+}$ when $A>-\delta$, and moreover, the BEBs at ${\mathbf{X}}_{BE-}$ can also be obtained according to the symmetry.

Definition 2.

If $\frac{\alpha (1+a)}{2\delta }>0$, i.e., $q<0$, then ${\mathit{X}}_{AE+}$ and ${\mathit{X}}_{PE}$, located at the opposite sides of ${\mathit{X}}_{BE+}$, may convert to each other via ${\mathit{X}}_{BE+}$, implying the occurrence of persistence; on the contrary, if $\frac{\alpha (1+a)}{2\delta }<0$, i.e., $q>0$, then ${\mathit{X}}_{AE+}$ and ${\mathit{X}}_{PE}$, located at the same side of ${\mathit{X}}_{BE+}$, may collide with each other at ${\mathit{X}}_{BE+}$ and disappear together, implying the occurrence of nonsmooth fold bifurcation.

3.3. Tangencies and the Visibilities

Tangencies and the related concepts are the essences of understanding the nonsmooth dynamics on discontinuity boundary, especially for sliding bifurcations of cycles [43,44]. Here we introduce the n-th order Lie derivatives ${\mathcal{L}}_{{\mathbf{F}}_{+}}^{n}H$ and ${\mathcal{L}}_{{\mathbf{F}}_{-}}^{n}H$ to discuss the tangencies as well as the local geometries on sliding boundaries $\partial {\mathrm{\Sigma }}_{S\pm }$(${\mathcal{L}}_{\cdots }^{n}H=\langle \frac{\partial {\mathcal{L}}_{\cdots }^{n-1}H}{\partial \mathbf{X}},\cdots \rangle$, where “⋯" denotes ${\mathbf{F}}_{+}$ or ${\mathbf{F}}_{-}$, $n\in N^{*}$).

Note that the two sliding boundaries in (8) are parallel to each other and two simple tangencies will exist on $\partial {\mathrm{\Sigma }}_{S\pm }$. Taking a point on the sliding boundary $\partial {\mathrm{\Sigma }}_{S+}$ as an example, there exist the following tangencies.

• Quadratic tangency, also be called the fold point, can be observed when the two conditions ${\mathcal{L}}_{{\mathbf{F}}_{+}}H=0$ and ${\mathcal{L}}_{{\mathbf{F}}_{+}}^{2}H\ne 0$ are satisfied, i.e.

  $$\left\{\begin{array}{c}\alpha y_{\mathrm{\Sigma }}+w+\delta =0,\hfill \\ -y_{\mathrm{\Sigma }}+z_{\mathrm{\Sigma }}\ne 0,\hfill \end{array}\right.$$

  (16)

  at which the trajectory driven by subsystem $S_{+}$ is tangent to sliding boundary $\partial {\mathrm{\Sigma }}_{S+}$. Furthermore, when ${\mathcal{L}}_{{\mathbf{F}}_{+}}^{2}H>0$, i.e., $-y_{\mathrm{\Sigma }}+z_{\mathrm{\Sigma }}>0$, it is a visible fold, while ${\mathcal{L}}_{{\mathbf{F}}_{+}}^{2}H<0$, i.e., $-y_{\mathrm{\Sigma }}+z_{\mathrm{\Sigma }}<0$, it is an invisible fold.
• Cubic tangency, also be called the cusp point, can be observed when the three conditions ${\mathcal{L}}_{{\mathbf{F}}_{+}}H=0$, ${\mathcal{L}}_{{\mathbf{F}}_{+}}^{2}H=0$ and ${\mathcal{L}}_{{\mathbf{F}}_{+}}^{3}H\ne 0$ are satisfied, i.e.

  $$\left\{\begin{array}{c}\alpha y_{\mathrm{\Sigma }}+w+\delta =0,\hfill \\ -y_{\mathrm{\Sigma }}+z_{\mathrm{\Sigma }}=0,\hfill \\ -\eta z_{\mathrm{\Sigma }}\ne 0,\hfill \end{array}\right.$$

  (17)

  at which one sliding orbit is tangent to sliding boundary $\partial {\mathrm{\Sigma }}_{S+}$. Furthermore, it is a visible cusp point when ${\mathcal{L}}_{{\mathbf{F}}_{+}}^{3}H·{\mathcal{L}}_{{\mathbf{F}}_{-}}H<0$, while it is a invisible cusp point when ${\mathcal{L}}_{{\mathbf{F}}_{+}}^{3}H·{\mathcal{L}}_{{\mathbf{F}}_{-}}H>0$.

Meanwhile, the tangencies and visibilities on $\partial {\mathrm{\Sigma }}_{S+}$ can be easily discussed according to the symmetry.

From the above analyses, we can find that a regular point ${\mathbf{X}}_{\mathrm{\Sigma }}$ on $\mathrm{\Sigma }$ can degenerate into a fold point via (16), and then a cusp point via (17), ultimately a BE via ${\mathcal{L}}_{{\mathbf{F}}_{+}}H={\mathcal{L}}_{{\mathbf{F}}_{+}}^{2}H={\mathcal{L}}_{{\mathbf{F}}_{+}}^{3}H=0$. On the other hand, we may also find that a cusp point separates the visibilities of the fold points while a BE can separate the visibilities of cusp points.

To date, the nonsmooth singularities on discontinuity $\mathrm{\Sigma }$ have been analyzed using nonconventional dynamics theory, which may help us to understand the dynamics mechanism for the nonsmooth structures observed both in the attractors of the fast subsystem and in whole system responses. In the following section, we push forward our research by numerical studies with certain parameter values. Furthermore, we discuss the dynamics generation mechanism of the non-bifurcation-induced transition structure characterized by the alternation from the outside spike attractor to the inner rest attractor in rest spike bistability.

4. Numerical Studies

When the values of the parameters other than the excitation amplitude A and the excitation frequency $\Omega$ are fixed, it is not difficult to understand that the excitation amplitude A determines the range of slow variable w while the excitation frequency $\Omega$, which has been assumed as $0<\Omega \ll 1$, does not. Based on the classical slow–fast decomposition method, the dynamics behaviors of the fast subsystem (4), such as attractors and their bifurcations, depend on the excitation amplitude A. Thus, when we slightly change the excitation frequency $Omega$, the response of the whole system may be constant. However, the situation is different when there is some special dynamical structure in (4).

In order to exhibit more details, we turn to numerical simulations via setting the parameters as the values in

Table 1. Parameters and their values.

Dimensionless Parameter	Value
$\alpha$	1.63
a	−0.3
b	0.5
$\delta$	−0.5
A	0.8
$\eta$	0.6

. From top to bottom in Figure 2, we give four different responses of the whole system (3) with different excitation frequencies.

When $\Omega =0.01$, only spiking oscillation mode can be observed, where the trajectory alternates between a cluster of relatively large amplitude oscillations (LAOs) and a cluster of relatively small amplitude oscillations (SAOs), repeated periodically. Then we can call this phenomenon mixed tonic spiking oscillations (MTSOs) since the amplitudes of spiking oscillations can be roughly divided into relatively large ones and small ones.

Interestingly when $\Omega =0.027$, the obvious QSs periodically appears between two SAOs, i.e., the so-called bursting phenomenon alternating between SPs and QSs can be observed in this case. Meanwhile, four alternations between LAOs and SAOs can also appear in every SP, leading to two clusters of LAOs and three clusters of SAOs in every SP. Then we can call such bursting phenomenon mixed bursting oscillations (MBOs).

Increasing the excitation frequency to $\Omega =0.027174$, MBOs still can be observed in this case. However, only two alternations between LAOs and SAOs can be observed in every SP, leading to that every SP consists of two clusters of SAOs separated by a cluster of LAOs.

Furthermore, by increasing the excitation frequency to $\Omega =0.0285$, alternations between SPs and QSs still exist in the oscillation structure. However, we may find that LAO no longer exists in SPs in this case, i.e., the amplitudes of SPs nearly keep constant. Therefore, we call this phenomenon typical bursting oscillations (TBOs).

Based on the numerical results in Figure 2, although the slight change of excitation frequency $\Omega$ does not affect the stabilities and bifurcations of attractors in the fast subsystem (4), the responses of whole system (3) from MTSOs to two different MBOs and then to TBOs appear in turn with increasing $\Omega$. We now try to reveal the generation mechanisms for those oscillation structures using the classical slow–fast decomposition method.

4.1. Generation Mechanism of MTSOs

Since the discontinuity boundary $\mathrm{\Sigma }$ is defined by $x=0$, we redraw the phase portrait in $(w,x)$ plane and the time history of x of the MTSOs when $\Omega =0.01$, respectively, given in Figure 3a,b. Then we may find that every spiking oscillation contacts the discontinuity boundary $\mathrm{\Sigma }$ (denoted as the blue horizontal line). Meanwhile, Figure 3a also reflects that the mixed tonic spiking oscillation with $\Omega =0.01$ may be not strictly periodic.

In order to exhibit more details about the oscillation structures, we choose ${\mathrm{\Pi }}_1:=w=0$, ${\mathrm{\Pi }}_2:=w=0.46$, ${\mathrm{\Pi }}_3:=w=0.55$ and ${\mathrm{\Pi }}_4:=w=0.78$ as cross sections, respectively, as shown in Figure 3a. Poincaré maps are plotted in Figure 4, where one may find that the response of the whole system (3) then is a quasi-periodic oscillation when $\Omega =0.01$, since all the Poincaré maps on the four cross sections are nearly closed curves.

According to the Poincaré maps in Figure 4, it is not difficult to predict that limit cycles, i.e., the so-called spike attractors exist in the Filippov fast subsystem (4). Furthermore, the relatively large closed curves intercepted from LAOs, labeled as LAO in Figure 4a–c, indicate that the corresponding spiking attractor behaves in crossing motions on $\mathrm{\Sigma }$. While the relatively small closed curves intercepted from SAOs, labeled as SAO in Figure 4b–d, indicate that the corresponding spiking attractor behaves in sliding motions on $\mathrm{\Sigma }$.

The relatively small closed curves intercepted from SAOs appears in the inner of the relatively large closed curves intercepted from LAOs, seen in Figure 4b,c. It indicates that the spiking attractor behaving in crossing motions coexists with the spiking attractor behaving in sliding motions. Moreover, the relatively small closed curves intercepted from SAOs behaves in different nonsmooth structures, which means that the corresponding spiking attractor may undergo sliding bifurcation.

Considering the symmetry, as shown by the green circles in Figure 3a, two spike attractors corresponding to the two SAOs exist in (4), respectively, labeled as $L{C}_{SA+}$ ($w>0$) and $L{C}_{SA-}$ ($w<0$). Taking $L{C}_{SA+}$ as the example, Figure 5 gives the corresponding geometry of sliding bifurcation, where $L{C}_{SA+}$ crosses through the visible cusp point $V{C}_{-}$ on $\partial {\mathrm{\Sigma }}_{S-}$ when $w\approx 0.533$ to undergo a multi-sliding bifurcation, seen in Figure 5b, leading to the single sliding segment in Figure 5a turns into two sliding segments in Figure 5c.

We now discuss the generation mechanism for the quasi-periodic MTSOs when $\Omega =0.01$ by using the classical slow–fast decomposition method. Via analyzing the conventional bifurcations of spike attractors in Figure 6a, we can give the complete details about bifurcations of spike attractors in

Table 2. Bifurcations of spike attractors as well as the abbreviations and parameter values.

Bifurcation	Abbreviation	Parameter Value
Fold bifurcation of limit cycles	$LP{C}_{2\pm }$	$w\approx \pm 0.566$
Multi-sliding bifurcation	$M{S}_{\pm }$	$w\approx \pm 0.533$
Fold bifurcation of limit cycles	$LP{C}_{1\pm }$	$w\approx \pm 0.444$

.

As shown by the red real arrows in Figure 6a, the relatively large amplitude spiking attractor, denoted as $L{C}_{LA}$, will collide with the two unstable cycles $L{C}_{US\pm }$, respectively, leading to the transitions from $L{C}_{LA}$ to $L{C}_{SA\pm }$ caused by the fold bifurcations of limit cycles $LP{C}_{2\pm }$. Similarly, the two relatively small amplitude spike attractors $L{C}_{SA\pm }$ may also collide with unstable cycles $L{C}_{US\pm }$, respectively, leading to the transitions from $L{C}_{SA\pm }$ to $L{C}_{LA}$ caused by the fold bifurcations of limit cycles $LP{C}_{1\pm }$.

Considering the symmetry, in order to exhibit the generation mechanism of the quasi-periodic MTSOs, we choose to superimpose one segment of trajectory with $\tau \in [\frac{3\pi +4n\pi }{2\Omega },\frac{5\pi +4n\pi }{2\Omega }]$ ($n\in N^{*}$) onto the one-parameter bifurcation diagram of spike attractors, seen in Figure 6b. During the slow variable w changes from $w=-0.8$ to $w=0.8$, it can be found that the evolution process of the trajectory can be explained well by using the bifurcation structures of spike attractors. For example, the multi-sliding bifurcation $M{S}_{-}$ leads the SAOs behaving in single-sliding oscillations to turn into a SAO behaving in double-sliding oscillations. Meanwhile, the two fold bifurcations of limit cycles $LP{C}_{1-}$ and $LP{C}_{2+}$ lead the transition from SAOs to LAOs and the transition from LAOs to SAOs, respectively.

It should be pointed out that the two multi-sliding bifurcations $M{S}_{\pm }$ lead to the change of sliding modes in SAOs rather than the transitions between attractors. We will ignore these two multi-sliding bifurcations because the main focus in our work is on the transition mechanism.

Since the whole system (3) is a slow–fast dynamical system when $0<\Omega <<1$, one may be interested in knowing whether or not the quasi-periodic MTSO response in Figure 3 should be a slow–fast dynamics behavior. The limit cycles in the fast subsystem (4) are presented in Figure 6a via giving the extreme points in the $(w,x)$ plane. Without loss of generality, we can also map those limit cycles as fixed points in the $(w,x)$ plane via discrete dynamics theory. As shown in Figure 6c, where $F{P}_{LA}$, $F{P}_{SA\pm }$ and $F{P}_{US\pm }$ denote the fixed points corresponding to $L{C}_{LA}$, $L{C}_{SA\pm }$ and $L{C}_{US\pm }$, respectively, while $S{N}_{1\pm }$ and $S{N}_{2\pm }$ denote saddle-node bifurcations corresponding to $LP{C}_{1\pm }$ and $LP{C}_{2\pm }$.

Similarly, we now represent the quasi-periodic MTSOs in Figure 3 with the maximum points of spikes, as shown by the cyan points in Figure 6c. Based on the classical slow–fast decomposition method, we can find that the quasi-periodic MTSOs is a relaxation oscillation induced by the four saddle-node bifurcations. Accordingly, we may say that the quasi-periodic MTSOs response in Figure 3 is a typical slow–fast dynamics behavior characterized by a relaxation oscillation on the torus.

4.2. Generation Mechanism of MBOs and TBOs

Based on the above analyses for the quasi-periodic MTSOs, the response of the whole system (3) can be well predicted by attractors and the bifurcation structures in (4). Meanwhile, according to the spike attractors and bifurcations in Figure 6a, it seems that the transitions between spike attractors should be the only alternative for the whole system trajectory. However, in the MBOs when $\Omega =0.027$, redrawn in Figure 7, the appearances of QSs connecting the neighboring two SAOs do not agree with the prediction based on the two fold bifurcations of limit cycles $LP{C}_{1\pm }$.

Figure 7 shows that the MBOs when $\Omega =0.027$ is a symmetrical period $3T_E$ response of system (3), where $T_E$ denotes the period of excitation $T_E=\frac{2\pi }{\Omega }$. For the convenience of analysis, we divide the whole period behavior into six equal time segments in time history according to the monotonous zones of slow variable w, labeled as $T_i(i=1,\dots ,6)$, respectively, seen in Figure 7b. It can be found that the QSs located on $\tau \in T_1$ and $\tau \in T_4$ represent the unique sticking motions which occur on discontinuity boundary $\mathrm{\Sigma }$. Since the sticking motion is an unique nonsmooth dynamics behavior that the trajectory in Filippov systems converges to a stable PE and behaves in a certain static state, we now count PE as well as BEBs into the bifurcation diagram in Figure 6a, as shown in Figure 8a.

By employing (9) and (15) as well as the results, it is easy to check that one stable focus type pseudo-equilibrium branch $PE$, existing in the sliding region ${\mathrm{\Sigma }}_S$, will disappear at the two boundary equilibrium points ${\mathbf{X}}_{BE\pm }$ via nonsmooth fold bifurcations (respectively labeled as $N{F}_{\pm }$). Under the separation of the four bifurcations $LP{C}_{1\pm }$ and $N{F}_{\pm }$, seen in Figure 8a, three special rest spike bistability structures, composed by stable focus type pseudo-equilibrium points (rest attractors) and nonsmooth stable cycles (spike attractors), can be observed in the fast subsystem. One is composed of $PE$ coexisting with the relatively large amplitude crossing cycle $L{C}_{LA}$, while the other two are composed of $PE$ and the small amplitude sliding cycles $L{C}_{SA\pm }$, more details are given in

Table 3. Rest spike bistability structures in the fast subsystem.

Sort	Attractor	Parameter Range	Border Bifurcation
I	$L{C}_{LA}$, $PE$	$w\in (-0.444,0.444)$	$LP{C}_{1-}$, $LP{C}_{1+}$
II	$L{C}_{SA+}$, $PE$	$w\in (0.444,0.533)$	$LP{C}_{1+}$, $N{F}_{+}$
III	$L{C}_{SA-}$, $PE$	$w\in (-0.533,-0.444)$	$LP{C}_{1-}$, $N{F}_{-}$

.

By then, the symmetrical period $3T_E$ MBOs in Figure 7 can be divided into two different transition routes. Taking the trajectories in decreasing intervals of the excitation as the example, as shown in Figure 8b, the segment with $\tau \in T_1$ (black dash–dotted curve) and the segment with $\tau \in T_3$ (black solid curve), respectively, converge to the special rest attractor $PE$ and the crossing spike attractor $L{C}_{LA}$. Then, the segment when $\tau \in T_1$ and the segment when $\tau \in T_3$ converge to the relatively small amplitude spiking attractor $L{C}_{SA-}$, respectively, caused by the bifurcation $N{F}_{-}$ and the bifurcation $LP{C}_{2-}$.

According to the bifurcation structures presented in the fast subsystem (seen in Figure 8a), even when starting from the stable PE, the trajectories of the whole system may also transition to $L{C}_{SA\pm }$ via the two nonsmooth fold bifurcations $N{F}_{\pm }$, and repeat the transition route similar with the first case with $\Omega =0.01$. Clearly, the transition route in the black dash–dotted curve is non-bifurcation-induced transition. Therefore, the key to understanding the generation mechanism of MBOs in Figure 7 is to reveal the dynamics mechanism underlying the appearance of the two coexisting distinct transition routes shown in the green circle in Figure 8b, especially the fundamental reason that results in the transition from $L{C}_{SA\pm }$ to $PE$.

One may note that the two coexisting attractors in the rest spike bistability II (III) may be very close to each other in the geometric position. For instance, see the rest spike bistability II when $w=0.445$ in Figure 9a. The clockwise stable cycle $L{C}_{SA+}$, starting from the smooth region $R_{+}$ ($R_{-}$), behaves in the sliding motion ${\mathrm{\Gamma }}_{S+}$ (${\mathrm{\Gamma }}_{S-}$) when it enters into ${\mathrm{\Sigma }}_S$ at the point $P_{R+}$ ($P_{R-}$), until it crosses through the sliding boundary $\partial {\mathrm{\Sigma }}_{S-}$ ($\partial {\mathrm{\Sigma }}_{S+}$) at the point $P_{S-}$ ($P_{S+}$) to enter into the opposite side of $\mathrm{\Sigma }$. We can find that the coexisting stable $PE$ is very close to the sliding segment ${\mathrm{\Gamma }}_{S-}$. We now focus on the division of attractor basins of the two coexisting attractors in the rest spike bistability II, which is important to explain the transition route from $L{C}_{SA\pm }$ to $PE$.

Since the stable $PE$ is very close to the sliding segment ${\mathrm{\Gamma }}_{S-}$ of spike attractor $L{C}_{SA+}$, the sliding trajectories (denoted by ${\Phi }_S(x\left(\tau \right),y\left(\tau \right),z\left(\tau \right))=(0,y_S\left(\tau \right),z_S\left(\tau \right))$) between $PE$ and the sliding segment ${\mathrm{\Gamma }}_{S-}$ will converge to one of the two coexisting attractors ultimately, According which the attractor basins of the two coexisting attractors in rest spike bistability II can be analyzed. Meanwhile considering that the sliding vector field (9) is linear, all the sliding trajectories ${\Phi }_S$ can be derived analytically according to that $PE$ is a stable focus, expressed by

$$\left\{\begin{array}{c}{y}_S\left(\tau \right)=\frac{{\mathrm{e}}^{\frac{-\tau }{2}}[{\mu }_1(5\sin \frac{\sqrt{35}\tau }{10}-\sqrt{35}\cos \frac{\sqrt{35}\tau }{10})+{\mu }_2(5\cos \frac{\sqrt{35}\tau }{10}+\sqrt{35}\sin \frac{\sqrt{35}\tau }{10})]}{6},\hfill \\ z_S\left(\tau \right)={\mathrm{e}}^{\frac{-\tau }{2}}({\mu }_1\sin \frac{\sqrt{35}\tau }{10}+{\mu }_2\cos \frac{\sqrt{35}\tau }{10}),\hfill \\ \frac{-w\pm \delta }{\alpha }\le y_S\left(\tau \right)\le \frac{-w\mp \delta }{\alpha },\hfill \end{array}\right.$$

(18)

where the two coefficients ${\mu }_{1,2}$ can be calculated by submitting the corresponding initial conditions of a sliding trajectory ${\Phi }_S$ into (18).

It is not difficult to verify that the two cusp points on sliding boundaries $\partial {\mathrm{\Sigma }}_{S\pm }$ both are visible, computed as $V{C}_{\pm }(0,-\frac{0.445\pm 0.5}{1.63},-\frac{0.445\pm 0.5}{1.63})$ via (17). Taking $V{C}_{+}$ as the initial conditions of (18), one special sliding orbit, labeled as ${\Phi }_{C+}$, can be computed with negative $\tau$. Certainly, the sliding trajectories between $PE$ and ${\Phi }_{C+}$ will slide to the stable $PE$ in the forward time. However, what is in doubt is whether all the sliding trajectories between ${\Phi }_{C+}$ and ${\mathrm{\Gamma }}_{S-}$ will converge to the spike attractor $L{C}_{SA+}$, that can be answered via the following numerical method.

One may note that all the sliding trajectories between ${\Phi }_{C+}$ and ${\mathrm{\Gamma }}_{S-}$ inescapably cross through the visible fold points between $V{C}_{+}$ and $P_{S+}$. Thus, the attractor basin structure can be easily obtained by which of the rest spike bistabilities II the trajectories from these visible fold points eventually converge to. Via numerical verification, the point $P_{B+}$ divides the sliding boundary $\partial {\mathrm{\Sigma }}_{S+}$ between $V{C}_{+}$ and $P_{S+}$ into two segments. As shown in Figure 9b,where trajectories starting from visible points between $V{C}_{+}$ and $P_{B+}$ converge to $PE$ to behave in sticking motions after a minimal oscillation in $R_{+}$, while trajectories starting from visible points between $P_{B+}$ and $P_{S+}$ oscillate to spike attractor $L{C}_{SA+}$. Similarly, another special sliding orbit, denoted as ${\Phi }_{sep}^{+}$, can be obtained numerically by taking $P_{B+}$ as the initial condition of (18). It can be regarded as an attractor basin boundary of the rest spike bistability II in sliding region.

Meanwhile, the classical slow–fast decomposition method reflects a fundamental fact in the family of slow–fast dynamical systems, referring to the fact that the dynamics behaviors of the fast subsystem are modulated continuously by the slow variable. According to that, the attractor basin structures of rest spike bistability II have changed after the whole system completes one spiking oscillation. Without loss of generality, taking the trajectory oscillating along with the $L{C}_{SA+}$ in Figure 9a as the example. Assume that the whole system trajectory starts at a point in the sliding region, where the corresponding slow variable is $W_1=0.5$sin$\left(0.027{\tau }_1\right)$, the trajectory completes a spike after one period of $L{C}_{SA+}$ (denoted as $T_{SA+}$). The slow variable has changed to be $W_2=0.5$sin$\left[0.027({\tau }_1+T_{SA+})\right]$, i.e., there exists a slight variation $\Delta w=|W_1-W_2|$. For more precision, here we name it the deviation of the slow variable (DSV for short).

Note that the excitation frequency has been assumed as $0<\Omega <<1$, then we have $\Delta w\approx 0.5\Omega T_{SA+}$, indicating that the DSV $\Delta w$ will increase along with the increase of $\Omega$, as shown in the sketch map in Figure 10a. As we have mentioned that $PE$ may be very close to the sliding segment ${\mathrm{\Gamma }}_{S-}$ of $L{C}_{SA+}$, seen in Figure 9a, the increase in DSV caused by the increase of $\Omega$, even though it is relatively small, can lead to unexpected influence on the attractor basins of rest spike bistability II. To illustrate this point, Figure 10b gives the attractor basin structures of rest spike bistability II with $w=0.445$, i.e., $\Delta w=0$. The special sliding orbit ${\Phi }_{sep}^{+}$ (colored in red) divides the sliding region ${\mathrm{\Sigma }}_S$ into two different sub-regions by the aid of sliding boundary ${\mathrm{\Sigma }}_{S+}$, respectively, denoted by ${\mathrm{\Sigma }}_{SP}$ and ${\mathrm{\Sigma }}_{SL}$. The sliding trajectories then converge to $PE$ when they are located in ${\mathrm{\Sigma }}_{SP}$ while the sliding trajectories converge to $L{C}_{SA+}$ (represented by the two sliding segments ${\mathrm{\Gamma }}_{S\pm }$) when they are located in ${\mathrm{\Sigma }}_{SL}$.

Interestingly, when we increase the DSV to a appropriate value, for instance, $\Delta w=0.034$, another unique sliding orbit besides ${\Phi }_{sep}^{+}$, denoted by ${\Phi }_{sep}^{-}$, can also be obtained numerically. Therefore, the attractor basin structures of rest spike bistability II can be redivided by using ${\Phi }_{sep}^{+}$ as well as the two sliding boundaries ${\mathrm{\Sigma }}_{S\pm }$, as shown in Figure 10c. Especially, we can find that the sliding segment ${\mathrm{\Gamma }}_{S-}$ then is located in ${\mathrm{\Sigma }}_{SP}$ while the sliding segment ${\mathrm{\Gamma }}_{S+}$ still is located in ${\mathrm{\Sigma }}_{SL}$. Sticking motions can be observed when the whole system trajectories oscillating along with $L{C}_{SA+}$ move to the sliding segment ${\mathrm{\Gamma }}_{S-}$ just in time. While spiking oscillations will persist when the whole system trajectories oscillating along with $L{C}_{SA+}$ move to the sliding segment ${\mathrm{\Gamma }}_{S+}$ just in time.

Accordingly, one may find that a threshold ${\Omega }_T$ exists in the whole system (3) under the parameter conditions in Table 1, which can be approximately computed as ${\Omega }_T\approx 0.025837$ via tracking the critical value of $\Omega$ of the appearance of the dash-dotted transition route in Figure 8b. Therefore, we have the following results for the transitions from rest spike bistability II to rest spike bistability I. When $\Omega <{\Omega }_T$, only the transition from $L{C}_{SA+}$ to $L{C}_{LA}$ can be observed, which is explained well by the fold bifurcation of limit cycles $LP{C}_{1+}$. However, when $\Omega \ge {\Omega }_T$, the special transition from $L{C}_{SA+}$ to $PE$ can also be observed besides the above. That can be explained by the disappearance of $L{C}_{SA+}$ induced by the nonnegligible DSV rather than a common bifurcation. Accordingly, we name it the transition induced by the DSV effect. Meanwhile, considering the symmetry, we have the same results about transitions from the rest spike bistability III to the rest spike bistability I.

For convenience, we can summarize the transition mechanism in a monotonous zone of slow variable w into the following two types when $\Omega \ge {\Omega }_T$.

• The regular transition route, meaning that the trajectory of the whole system will transition to $L{C}_{LA}$ from $L{C}_{SA\pm }$, and then to $L{C}_{SA\mp }$ via the fold bifurcations of limit cycles, leading to the birth of mixed tonic spiking oscillation pattern;
• The irregular transition route, meaning that the trajectory of the whole system will transition to $PE$ from $L{C}_{SA\pm }$ via DSV effect, and then to $L{C}_{SA\mp }$ via the nonsmooth fold bifurcations, leading to the birth of bursting oscillation pattern.

By then, the generation mechanism underlying the symmetrical period $3T_E$ MBOs in Figure 7 can be understood well based on the above analyses. More precisely, assuming that the trajectory starts at the maximum of w in $T_1$, it, respectively, passes through the rest spike bistability II and III three times when it oscillates along with spike attractor $L{C}_{SA\pm }$, where the transition induced by DSV effect, respectively, only takes place once. It results in a symmetrically periodic transition structure composed of two irregular transitions followed by for regular transitions in a complete oscillation period. Then, the period $3T_E$ MBOs, consisting of two QSs characterized by sticking motions and two SPs characterized by mixed spiking oscillations, can be observed in the whole system.

Based on the analysis of the generation mechanism of the period $3T_E$ MBOs when $\Omega =0.027$, it is not difficult to understand the responses of whole system (3) when $\Omega =0.027174$ and $\Omega =0.0285$. We now focus on the generation mechanisms of the two responses in Figure 11.

When $\Omega =0.027174$, seen in Figure 11a, we may find that the whole system (3) behaves in an asymmetrical period–$1T_E$ MBOs, and the trajectory, respectively, passes through the rest spike bistability II and III once when it oscillates along with spike attractor $L{C}_{SA\pm }$. The transition induced by the DSV effect only occurs when the trajectory moves to the rest spike bistability I from the rest spike bistability III, however, it does not take effect when the trajectory moves to the rest spike bistability I from the rest spike bistability II. It results in the appearance of the asymmetrical period–$1T_E$ MBOs consisting of one QS characterized by a sticking motion and one SP characterized by mixed spiking oscillations.

When $\Omega =0.0285$, seen in Figure 11b, the whole system response is a symmetrical period–$1T_E$ TBOs. Similarly, with the whole system response with $\Omega =0.027174$, the trajectory in this case also passes through the rest spike bistability II and III once when it oscillates along with spike attractor $L{C}_{SA\pm }$, respectively. However, we can find that the transition induced by the DSV effect may take effect at the two transition locations, leading to the trajectory transition from $L{C}_{SA\pm }$ to $PE$. That implies the appearance of the symmetrical period–$1T_E$ TBOs consisting of two QSs characterized by sticking motions and two SPs characterized by spiking oscillations without LAO.

4.3. Coexistence of Whole System Responses

According to the generation mechanism of whole system responses involving transitions induced by the DSV effect, one may find that any response of the whole system is a certain combination form of regular and irregular transition routes. Moreover, such a combination may be random in theory since period–3 movement has been obtained when $\Omega =0.027$, implying the appearance of an arbitrary period of whole system response, even chaos.

More interestingly, the period–$1T_E$ MBOs in Figure 11a is asymmetrical, i.e., there must be another asymmetrical period–$1T_E$ MBOs coexisting with the one in Figure 11a based on the symmetry. It indicates that the coexistence of attractors can also be observed in the whole system. We will give further discussions in this section.

Figure 12a gives the coexistence response of the whole system when $\Omega =0.027174$, and the corresponding initial conditions from top to bottom are $\left(x\right(0),y(0),z(0\left)\right)=(0,0.5,0.5)$, $\left(x\right(0),y(0),z(0\left)\right)=(0,-0.5,-0.5)$ and $\left(x\right(0),y(0),z(0\left)\right)=(-2,-2,2.2)$, respectively. Numerical results show that there also exists a symmetrical period–$1T_E$ MTSOs (seen in Figure 12(a1)) besides the two coexisting asymmetrical period–$1T_E$ MBOs (seen in Figure 12(a2,a3)), indicating that tristability then exists in the whole system (3). Noting that the asymmetrical period–$1T_E$ MBOs in Figure 12(a3) is symmetrical with the one in Figure 11a, the corresponding generation mechanism can be obtained according to the symmetry. Meanwhile, for the symmetrical period–$1T_E$ MTSOs in Figure 12(a1), the generation mechanism is attributed to the fact that only two regular transition routes appear in a complete period motion.

In order to show more details about the coexistence of responses of the whole system (3), we turn to exhibit the responses of the whole system via plotting the Poincaré map. The whole system trajectory will oscillate with $L{C}_{LA}$ when it is in the regular transition route. In contrast, the whole system trajectory may move along with $PE$ in the irregular transition route. Accordingly, ${\mathrm{\Pi }}_1:=w=0$ is chosen as the Poincaré cross section, which can reflect the symmetry of the whole system response and can distinguish well between the two transition routes. As shown in Figure 12b, we give the Poincaré map of the whole system when $\Omega \in (0.02685,0.027183)$, in which the period–$1T_E$ MTSOs (colored in orange) can persistently exist in the whole frequency range.

Meanwhile, the period–$1T_E$, period–$3T_E$, period–$2T_E$, period–$5T_E$ and period–$3T_E$ MBOs can discontinuously appear in turn along with the decrease of $\Omega$. In those periodic MBOs, further numerical simulations show that the irregular transition route will appear at most once, which is absent here for brevity. Here the periodic MBOs colored in red denote the irregular transition route that only occurs nearby the rest spike bistability II. While the ones colored in green denote the irregular transition route that only occurs nearby the rest spike bistability III. Moreover, the black ones denote the irregular transition route occurs both in the rest spike bistability II and in III. Furthermore, the structures of MBOs may also alternate between asymmetry and symmetry in this progress. It leads the coexistence of responses of the whole system (3) to alternate between tristability and bistability, more details see

Table 4. Coexistence of the responses of whole system (3).

Coexistence Type	$\mathbf{\Omega }$ Zone	Response Type
tristability	$(0.0271193,0.027183)$	one period–$1T_E$ MTSOs 1 and two period–$1T_E$ MBOs 2
bistability	$(0.0269777,0.0270228)$	one period–$1T_E$ MTSOs 1 and one period–$3T_E$ MBOs 1
tristability	$(0.0269256,0.0269444)$	one period–$1T_E$ MTSOs 1 and two period–$2T_E$ MBOs 2
bistability	$(0.0269029,0.0269114)$	one period–$1T_E$ MTSOs 1 and one period–$5T_E$ MBOs 1
tristability	$(0.0268919,0.0268956)$	one period–$1T_E$ MTSOs 1 and two period–$3T_E$ MBOs 2

¹ Symmetrical response(s). ² Asymmetrical response(s).

.

Then, we can find that the whole system (3) exhibits various dynamics behaviors characterized by one period–$1T_E$ MTSOs coexisting with one (or two) periodic MBOs when the irregular transition route is involved. It indicates that the response only involving regular transition route still exists in the whole system when $\Omega \ge {\Omega }_T$. We now turn to show another interesting coexistence of the whole system responses.

When $\Omega =0.0285$, tristability structure can also be observed in the whole system. As shown in Figure 13a, one symmetrical period–$1T_E$ TBOs, which has been presented in Figure 11b, can coexist with two asymmetrical period–$1T_E$ MTSOs, the corresponding initial conditions from top to bottom are $\left(x\right(0),y(0),z(0\left)\right)=(0,0,0)$, $\left(x\right(0),y(0),z(0\left)\right)=(2,2,-2.2)$ and $\left(x\right(0),y(0),z(0\left)\right)=(-2,-2,2.2)$, respectively. Similarly, we also give the Poincaré map when $\Omega \in (0.0284,0.02852)$ in Figure 13b, where the symmetrical period–$1T_E$ TBOs (colored in orange) can persistently exist in the whole frequency range. Meanwhile, both the two asymmetrical period–$1T_E$ MTSOs presented in Figure 13(a2,a3) perform the obvious period-doubling cascades along with the decrease of $\Omega$, colored in light gray and gray, respectively.

Different from the former case in Figure 12b, one may find that tristability structure always exists in the whole system (3) when $\Omega \in (0.0284,0.02852)$. Moreover, the tristability can be composed not only by a period–$1T_E$ TBOs coexisting with two periodic MTSOs but also by a period–$1T_E$ TBOs coexisting two chaotic MTSOs.

5. Conclusions and Discussions

Based on the unique rest spike bistability consisting of a limit cycle and the inner pseudo equilibrium in a Filippov slow–fast dynamical system, we have further discussed the non-bifurcation-induced transition characterized by the transition from the outside spike attractor to the inner rest attractor via various slow–fast dynamic behaviors in system (3). We found that the common feature DSV effect induced by the objective reality of timescale difference is one key factor in the transition. Furthermore, it has been proved that the DSV effect can be determined by timescale difference. In this sense, the non-bifurcation-induced transition mechanism can be regarded as the comprehensive effect of rest spike bistability and timescale difference.

Our work also reveals that there exists an essential threshold in timescale difference when a slow–fast dynamical system possesses such an interesting transition mechanism, which can be used to distinguish between transition routes across the whole system responses. As in the Filippov slow–fast dynamical system (3), the threshold ${\Omega }_T$ under the parameter conditions in Table 1 has been given numerically. The whole system only behaves in a regular transition route that can be predicted well by bifurcation structures when $\Omega <{\Omega }_T$. However, the irregular transition route, i.e., the non-bifurcation-induced transition route may coexist with the regular transition route when $\Omega >{\Omega }_T$. The whole system responses then tend to be complicated and varied (shown in Figure 12b and Figure 13b). This manifests mainly in the two aspects as follows.

• The combination form of the two coexisting transition routes may be random since the period-3 solution can be observed in the whole system responses, leading to various oscillation patterns such as multi-periodic and chaotic slow–fast dynamics behaviors.
• Interesting and various coexistences of whole system responses can also be observed. For instance, discontinuous alternations between bistability, tristability, and persistent tristability may be performed across the whole system.

Note that the core of the classical slow–fast decomposition analysis method extensively adopted in slow–fast dynamics represents the dynamics analysis of the fast subsystem after timescale separation via the key slow variable as the bifurcation parameter. The influence of timescale difference on the whole system is ignored. On the contrary, our results show that the timescale difference must be considered in certain conditions, such as the appearance of rest spike bistability in this paper. Moreover, our numerical results indicate that the timescale difference is highly sensitive to the whole system in this regard, and thus warrants the attention of researchers in fields related to slow–fast dynamics.

A similar non-bifurcation-induced transition as that used in this paper has been involved in recent works [22,31]. However, to the best of our knowledge, we are the first to investigate that such a special transition mechanism should be a comprehensive effect of rest spike bistability and timescale difference. Furthermore, our work suggests that one non-bifurcation parameter of the fast subsystems, i.e., the excitation frequency, can also induce various slow–fast dynamic behaviors in a Filippov slow–fast dynamical system. It may provide a valuable reference for the modification and design of artificial neuron models since nearly all existing artificial neuron models are discontinuous.

Additionally, this paper mainly introduces the non-bifurcation-induced transition mechanism under the comprehensive effect of rest spike bistability and timescale difference via numerical methods. An in-depth investigation into the innate character underlying this transition mechanism and to theorize the related numerical results will prove meaningful for understanding slow–fast dynamics, which is not a trivial problem and will be our primary concern in the future.