Multi-bifurcation cascaded bursting oscillations and mechanism in a novel 3D non-autonomous circuit system with parametric and external excitation

In this paper, multi-timescale dynamics and the formation mechanism of a 3D non-autonomous system with two slowly varying periodic excitations are systematically investigated. Interestingly, the system shows novel multi-bifurcation cascaded bursting oscillations (MBCBOs) when the frequency of the two excitations is much lower than the mean frequency of the original system. For instance, periodic, quasi-periodic and chaotic bursting oscillations induced by a variety of cascaded bifurcations are first observed, and the phenomenon of spiking transfer is also revealed. Besides, stability and local bifurcations of the system are comprehensively investigated to analyze the mechanism of the observed MBCBOs, in which bifurcation diagram, Lyapunov exponents, time series, phase portraits, and transformed phase diagrams are used. Finally, through a circuit simulation and hardware experiment, these complex dynamics phenomena are verified physically.


Introduction
Multi-timescale dynamics are ubiquitous in diverse research fields. For example, by regulating calcium and potassium activity, the biophysical model exhibits the spontaneous bursting behavior in the developing retina [1], complex mixed-mode oscillatory patterns are revealed in a Belousov-Zhabtinsky reaction with fastslow timescale [2], the coexisting fast-scale and slowscale bifurcation are observed in simple dc/dc converters [3] and the complicated bursting oscillations are generated in the hydroelectric power system with multiple time scales [4]. In 1963, Lorenz discovered the first chaotic system when studying the weather model [5]. Thereafter, various chaotic systems have been proposed by many scholars, especially Lorenz system family. As a type of classically chaotic system, Lorenz system family has been extensively studied and applied in the past decades. However, most of the work focused on the study of dynamic characteristics at the same timescale [6][7][8], and the multi-timescale dynamics of Lorenz system family is rarely reported. Therefore, "what kind of dynamic behavior can Lorenz system family with multi timescale and multi-frequency periodic excitation exhibit?" remains an open question that is worth studying further.
Bursting oscillations (BOs) caused by periodic switching of the spiking state (SP) and the quiescent state (QS) [9], as an interesting multi-timescale dynamics, has great potential for applications. For example, researchers in [10] used several types of BOs to drive electrodynamical systems. Authors in [11] demonstrated a specific type of artificial neuron with threshold-driven spiking and built a neural network for pattern recognition. Furthermore, it can also be applied to the study of stability and control design of pulsed power supply systems according to similar oscillation characteristics [12]. Therefore, the exploration of new patterns of BO can not only supplement the existing theory but also provide a variety of options for its application. Over the past few decades, a large number of BOs with different bifurcation mechanisms have been investigated. For instance, in some two-timescale system [13][14][15][16], fold/fold, fold/Hopf and symmetric Hopf/Hopf burster have been observed. Recently, more complex BOs are found on multi-timescale systems. Such as the turnover-of-pitchfork-hysteresis-induced and the compound pitchfork-hysteresis-induced bursting oscillations are demonstrated in [17]. Otherwise, some complex mixed-mode oscillations (MMOs) are revealed in [18,19]. The discovery of this kind of bursting oscillation provides a new direction for research. Although a substantial amount of research works has been carried out for multi-timescale systems, many complex BOs or MBCBOs still remain unexplored. This inspires our work on the analysis of dynamics of multi-timescale.
Based on the above discussion, this work is motivated to investigate the multi-timescale dynamics of a novel parametrically and externally driven circuit system, which the original system belongs to the class of the Lorenz system family. It is noteworthy that the system exhibits complex MBCBOs. For instance, periodic, quasi-periodic and chaotic MBCBOs are uncovered, which are caused by the cascading of pitchfork bifurcation-delay-induced hysteresis (PBDIH) [20], homoclinic orbit (HO) and Hopf bifurcation (HB). Unfortunately, these complicated MBCBOs were never studied, and the mechanism of various modes of MBC-BOs is not revealed in detail. Noticeably, the period of a single MBCBOs is determined by excitation at a small-frequency. Moreover, the effect of system parameters on dynamic is also thoroughly studied by adjusting the system parameters. Finally, an analog circuit is Table 1 If f 1 < f 2 , the relation between g(t) and f (g(t))

Table 2
If f 1 > f 2 , the relation between g(t) and f (g(t)) designed to verify the above complex oscillation signal through circuit simulation and hardware experiments. The rest of this paper is organized as follows. In Sect. 2, a non-autonomous system with parametric and external excitation is presented and the excitation transformation method is briefly introduced. In Sect. 3, stability and local bifurcation of the proposed 3D nonautonomous system are comprehensively investigated. In Sect. 4, the complex dynamic behaviors are studied by numerical simulation. Besides, circuit simulations and hardware experiments are performed to validate the numerical simulation in Sect. 5. Finally, conclusions are drawn in Sect. 6.

Mathematical model
In 2007, Cai et al. obtained a new Lorenz-like chaotic system by modifying the third equation of Lorenz system [21], studied its basic dynamic characteristics in detail, and designed a single controller to control the system. In fact, this system inherits basic characteristics of Lorenz system, such as large Lyapunov exponent and butterfly chaotic attractor. Therefore, we choose the system as research object to study the multi-timescale dynamics of Lorenz system family, and its model is described by where x, y, and z represent three state variables, while a, b, and c are positive constants. Based on system (1), a novel non-autonomous system with parametric and external excitation is proposed to model multitimescale dynamics of Lorenz system family under periodic excitation. Different from the previous systems, in this system, parametric excitation and external excitation are introduced separately to two differential equations. Its mathematical model is depicted as, where w 1 = β 1 cos(2π f 1 t) and w 2 = β 2 cos(2π f 2 t) are two excitation signals. If the variable of system (1) is regarded as a harmonic signal with continuous phase and random envelope change, then the mean frequency of the original system defines the mean time de rivative of the phase [22,23], as follows where N is the number of maxima or minima in the time T . Here, the frequency ratio between fast subsystem and slow subsystem is defined as the order gap σ =f o / f s , wheref o and f s is the mean frequency of the original system and the excitation frequency, respectively. When the system (1) parameters are set as a = 10, b = 2 and c = 1, theω o can be calculated as 16.71, which means thatf o =ω o /(2π) ≈ 2.659. As a result, when the frequency of the two excitation in the system (2) meets 0 < f s 1 (that is, far less thanf o ), the order gap of the improved non-autonomous system is much larger than 1, which inevitably leads to the generation of multi-time scale dynamics.

Methodology
For a fast-slow dynamical system with two slowly varying excitation, its dynamic behavior can be described bẏ where vector x(t) ∈ R n stands for fast subsystem. Usually, w 1 (t) = β 1 cos(2π f 1 t) and w 2 (t) = β 2 cos(2π f 2 t) is regarded as a slow subsystem. It is considered that the classical fast-slow analysis method [24] cannot be used to analyze the dynamics of the fast subsystem when there are multiple variables in the slow subsystem. Therefore, to address this problem, it is necessary to introduce a generalized slow variable g(t) by the excitation transformation method [18], so that w 1 (t), w 2 (t) and g(t) exist uniformly in the system in the form of independent variables and dependent variables of the function. Next, brief introduction of the excitation transformation method is given as below. Assuming that two slow variables are w 1 = cos(at) and w 2 = cos(bt), and n = a/b, where n is integer, rational or irrational number. Let's take n as an integer to illustrate the excitation transformation in detail. For the other two cases, [18,25] give the detailed explanations. Firstly, according to Moivre's formula [26], the following relationships can be obtained, Then, by expanding the left side of polynomial (5), we can get Comparing the right side of (5) and (6), we have Using trigonometric Eq. (7) can be rewritten as, where m ≤ n and m is even.
If there exists an integer k, such that a = pk and b = qk are true. Depending on (8), we can get It follows that the periodic excitation w 1 and w 2 in the system can be replaced by a polynomial function of cos(kt). It turns out that there is only a uniformly generalized slow variable in the system, cos(kt), denoted by g(t).
Secondly, according to different excitation frequencies, the system can be divided into slow subsystem g(t), and generalized fast subsystem in the followling three cases: , then fast subsystem can be rewritten as , then fast subsystem can be rewritten as It is concluded that slow subsystem can be replaced by g(t) with only one slow-varying excitation. In addition, the excitation amplitude β 1 and β 2 are usually taken as a fixed parameters when performing fast-slow analysis. Note that here the excitation amplitude is generally 1, and typical examples are given in Tables 1 and  2. Consequently, it can be seen from (9) that the highfrequency excitation can be replaced by the polynomial of low-frequency excitation. Next, these complicated dynamics are analyzed theoretically and numerically.

Local bifurcation and stability analysis of equilibrium point
For a dynamic system, the characteristics of the equilibrium point affect the complex dynamics to a large extent. For example, unstable saddle-focus may lead to the formation of chaotic attractors [27], the changes of equilibrium branch structure and properties lead to the emergence of all kinds of bifurcations [28]. Therefore, this section will focus on the characteristics of the equilibrium point and the local bifurcation. First of all, suppose thaṫ it is obvious that the equilibrium points of system (2) are determined by the following formula where q = cw 1 + w 2 + bc. Therefore, it is necessary to analyze the existence of equilibrium points in two cases: It can be found that the existence of the equilibrium point is closely related to w 1 and w 2 , which means that the trajectory of equilibrium point (TEP) will undergo a complicated change when the two slowly varying excitations change at the same time.
The local stability of the equilibrium point can be determined by eigenvalues of following Jacobian matrix The characteristic equation of (10) can be found as According to the Routh-Hurwitz criterion [29], if the following conditions are met Otherwise, it is unstable. Now, a specific example is given to illustrate the characteristics of equilibrium change. Letting a = 10, b = 2, c = 1, β 1 = β 2 = 3, f 1 = 0.003 and f 2 = 0.001, then w 2 = g(t), w 1 can be replaced by 4g(t) 3 − 3g(t) based on Table 2. Next, the corresponding equilibrium points E 0 and E ± are substituted into (11), and the following characteristic equation can be obtained as, and According to the characteristic equations of (16) and (17) with higher degree polynomials, it can be known that there are many local bifurcations under the control of the slowly varying excitation w 1 and w 2 . And there are mainly two kinds of bifurcations: (i) Pitchfork bifurcation [30].
On the one hand, if the following condition is true, i.e. g(t) = −0.836, the characteristic equation has a zero eigenvalue, which proves that this equilibrium point is a non-hyperbolic equilibrium point. When g(t) passes through -0.8362 slowly, the stable equilibrium (SE) point is split into two SE points and one unstable equilibrium (UE) point, leading to the emergence of pitchfork bifurcation (BP). The corresponding TEP is plotted in Fig. 1c, it also can be seen that the bifurcation point is supercritical.
On the other hand, if the following condition is satisfied, i.e. g(t) = 0.815, −0.455 and −0.360, the characteristic equation has a pair of purely imaginary eigenvalues, which implies that the TEPs and periodic solutions meet at a bifurcation point. Thus, Hopf bifurcation may occur, which leads to periodic oscillation and unstable limitcycle (UL). Besides, its period is determined by the imaginary part of the eigenvalue. Following the above analysis, the TEPs under the condition of other parameters are also plotted in Figs. 1 and 2 by utilizing XPPAUT. It is verified that the system has a variety of TEPs for different system parameters and the excitation frequency. Besides, the local bifurcation points have an increasing trend with the increase of frequency ratio. This behavior is likely to lead to complex dynamical behavior.

Numerical results of dynamics and fast-slow analysis
In this section, a variety of complex oscillatory behaviors of system (2) will be simulated by using MATLAB. Besides, the ODE45 algorithm is used to solve differential equations, in which the time step is set as 0.01. Furthermore, various MBCBOs of system (2) under the different frequency ratio are studied separately. And the effect of system parameters on dynamical behaviors is also studied.

The delayed-pitchfork/Hopf/homoclinic induced periodic bursting
Setting f 1 = f 2 = 0.001, β 1 = β 2 = 3, a = 10, b = 2, c = 1, and the initial condition is (0.1, 0, 0). The simulation time is set as 5000. Under this parameter condition, the generalized fast subsystem is Eq. (12), and the generalized slow subsystem is g(t). Dynamic and fast-slow analyses are presented in Fig. 3. It can be seen from Fig. 3b that the evolution of BO and equilibrium point is basically consistent, and the intersection of Hopf bifurcation set (HBS) and E ± represents HB point.
The formation mechanism of BO is further analyzed by combining Fig. 3b, c. Suppose the trajectory starts at QS 1 (i.e. the point A in Fig. 3c) under the control of generalized slow variable g(t). When the trajectory passes through BP slowly, the pitchfork bifurcation-delay-induced hysteresis is formed, which results in the trajectory going through a period delay before entering SP 1 . Then, because the equilibrium point is always stable, the trajectory gradually converges to QS 2 . Following the trajectory passes through the HB + point slowly, which leads to SP 2 and the amplitude increases gradually. At the same time, two symmetrical saddle-focus appear in the system. Therefore, HO is formed, and the trajectory begins to oscillate back and forth between the two saddle-focus. Finally, the trajectory returns to HB + again, leading the trajectory returns to QS 1 until BP. Consequently, the oscillation behavior is named as the delayed-pitchfork/Hopf/homoclinic induced periodic bursting.

The first cascading periodic bursting induced by delayed-pitchfork/Hopf/homoclinic
Leave the other parameters unchanged and increase f 1 to 0.002. Under this parameter condition, the generalized fast subsystem is Eq. (10), and the generalized slow subsystem is g(t). Then dynamic and fast-slow analysis is plotted in Fig. 4. It is seen that the period of a single cascaded bursting oscillation is determined by excitation at a small-frequency and is equal to T = 1/ f 1,2 . Figure 4 shows that there is a cascade of two patterns within a single bursting oscillation period. In fact, there is one more PBDIH induced bursting [17] comparing with the delayed-pitchfork/Hopf/homoclinic induced periodic bursting. Therefore, the mechanism is similar to the former, refer to the previous subsection.

The second cascading periodic bursting induced by delayed-pitchfork/Hopf/homoclinic
Increase continuously f 1 to 0.003. Dynamic and fastslow analysis is shown in Fig. 5. Besides, Fig. 1c shown that the TEP has changed tremendously comparing with the previous two patterns. The most prominent feature is the bending of the TEP, and there are three pairs of HB. Actually, the mechanism of bifurcation is the same as The first cascading periodic bursting induced by delayed-pitchfork/Hopf/homoclinic. From Fig. 5b, c we can also know that when the trajectory gets to the BP point, it experiences a PBDIH before entering SP. The main difference is that there are continuous twists and turns, and the state of system (2) does not change at all when the trajectory passes through the HB1 + and HB2 + points, and stay in QS. Subsequently, system (2) enters the periodic SP because the equilibrium is unstable. Finally, trajectory passes through the HB3 + , the trajectory returns to QS. Therefore, the oscillation is also considered as the second cascading periodic bursting induced by delayed-pitchfork/Hopf/homoclinic.

The first cascading quasi-periodic bursting induced by delayed-pitchfork/Hopf/homoclinic
A large number of numerical results show that the effects of system parameters on the system dynamics are very similar under different frequency ratios. In this subsection, two determined frequencies are selected to study the system dynamics. Now setting a = 10, c = 1, β 1 = β 2 = 3, f 1 = 0.002, f 2 = 0.01, and the initial condition is (0.1, 0, 0). The simulation time is set as 5000. Under this parameter condition, the generalized fast subsystem is Eq. (11), and the generalized slow subsystem is g(t). Besides, w 1 = g(t), w 2 = f (g(t)) = 16g(t) 5 − 20g(t) 3 + 5g(t). Then the two-parameter local bifurcation sets for g(t) and b can be obtained, and shown in Fig. 6. It is seen that the number and location of BP and HB will change with b and g(t).

Bifurcation diagram and Lyapunov exponents analysis
To further study the effect of parameter b, other parameters are fixed and b is used as a variable parameter to draw the bifurcation diagram and Lyapunov exponents of single parameter, results are shown in Fig. 7. It can be seen that the dynamical behaviors reflected by the trajectories in Fig. 7a, b agree with each other, and system (2) exhibits complex oscillatory behaviors.  3). Subsequently, in the range of (3, 4.4), (8.3, 9) and (4.4, 8.3), the system exhibits quasi-periodic and chaotic oscillation respectively, which means that the system leads to chaos through the quasi-periodic route. Finally, the system exits chaos by inverse quasi-periodic. Furthermore, it is obvious that the trajectory jumps around the zero-point within the region (−4, 8), which means that complex bursting oscillations may occur. Setting b = 2 and using g(t) as a variable parameter, the bifurcation diagram and Lyapunov exponents of single parameter are shown in Fig. 8. The results show that under the control of the slowly varying parameter g(t), the system experiences all kinds of attractors, including periodic, quasi-periodic and chaotic attractors. Besides, some interesting bifurcations have also been observed, such as BP, HB, forward period-doubling bifurcation (FPDB) and reverse period-doubling bifurcation (RPDB). Therefore, it is also confirmed that the bursting oscillation is caused by the switching between different attractors. Also, it is interesting to observe from LEs in Fig. 7b that the overall behavior of the system presents a quasi-periodic state. However, from the time series in Fig. 9a, it seems to be in chaos.

Mechanism analysis
For a better understanding of the bursting process, Dynamic and fast-slow analysis are presented in Fig. 9. It is shown from Fig. 9b, c that complex phenomenon of spiking transfer under the control of HB. Moreover, the trajectory of purple, light blue, red and brown represents four different patterns of spiking in Fig. 9d.
Assuming that the trajectory starts to oscillate at point A. Firstly, the trajectory experiences the first pattern of BO caused by PBDIH when it passes through BP 1 and the oscillation trajectory is represented by the purple. However, the second pattern of BO is more complicated comparing with the former. When g(t) moves to BP 3 , the trajectory enters the second pattern of BO, at this time, the SP 2 is formed. The existence of HB 1 causes the saddle-foci equilibrium point to appear, and then the trajectory oscillates between two saddle-foci and is connected by HO. Finally, the system returns to QS through HB 2 . Next, the trajectory follows the SE passes through HB 3 , and the trajectory enters into SP 3 . Due to the existence of HO, the trajectory oscillates between two saddle-focus again. Subsequently, the trajectory returns to QS by HB 3 again. Then, the pattern spiking SP 3 is genetared by HB 2 . Although both are caused by HB, the amplitude of oscillation of the third model is significantly smaller than that of the first two, and oscillates in the lower half of the equilibrium point. Sequentially, the system exits SP 3 return to QS at HB 1 . Finally, the trajectory follows the equilibrium point to BP 2 and returns to the first pattern of BO. So far, the whole oscillation process has been analyzed completely. This type is also called delayed-pitchfork/Hopf/homoclinic cascade induced quasi-periodic bursting .because it keeps switching between these bifurcations.

The second cascading quasi-periodic bursting induced by delayed-pitchfork/Hopf/homoclinic
To further explore the effect of parameter b, we choose b = 4 to analyze the system dynamic.

Bifurcation diagram and Lyapunov exponents analysis
The bifurcation diagram and Lyapunov exponents with varying g(t) can be acquired from system (2), and shown in Fig. 10. Interestingly, the local bifurcation mainly occurs in the region of g(t) < 0, in which the BP and three pairs HB are located at

Mechanism analysis
The time series of x and g(t) are presented in Fig. 11a, b. While the superposition graph for the evolution of equilibrium point and time series of x and transformed phase diagram are plotted in Fig. 11c, d. It is obvious that the oscillation process can be divided into three parts, which are represented by green, pink and purple trajectory respectively. Firstly, the pattern 1 is induced by cascaded PBDIH-HO-HB. Unlike the others, the pattern 2 is induced by HB-HO-HB, and the pattern 3 is induced by cascaded HB-HO-HB-PBDIH. Therefore, such oscillation is also known as a quasi-periodic bursting caused by PBDIH-HO-HB/HB-HO-HB/HB-HO-HB-PBDIH cascade.

4.6
The cascading chaotic bursting induced by delayed-pitchfork/Hopf/homoclinic When b increases to 5, a novel class of chaotic bursting can be observed from system (1).

Bifurcation diagram and Lyapunov exponents analysis
The bifurcation diagram and Lyapunov exponents of generalized slow variable g(t) are shown in Fig. 12. From which, we can konw that the path for system (1) to enter the chaotic state is HB, FPDB. The system exhibits chaos in most regions with slow-varying parameter. Besides, several period windows are generated by system (2). At the same time, Fig. 7b shows that the maximum Lyapunov exponent LE1 = 0.2154, which means that the motion is chaotic.

Mechanism analysis
The time series of x are plotted in Fig. 13a, b. There are two kinds of spiking transitions in a SP, including periodic spiking states (PSP) in the region II and IV, and chaotic spiking states (CSP) in the region I, III and V. The superposition graph for the evolution of equilibrium point and time series of x, and transformed phase diagram are presented in Fig. 13c, d separately. Firstly, when the trajectory slowly passes through the BP point, the system (1) enters the CSP through PBDIH. Then, the trajectory reaches the HB2 ± point, which generates PSP from the upper and lower branches of the equilibrium. Next, the trajectory exits the PSP and entering CSP again by HB3 ± point. Finally, the trajectory returns to QS by BP point again. Therefore, the oscillation process can be understood as a chaotic bursting caused by delayed-pitchfork/Hopf/homoclinic cascade.

Circuit design and hardware experiment
In this section, the above dynamics are realized and verified by circuit simulation and hardware experiment. The simulation is completed by the Multisim 13.0 circuit simulation platform, and the hardware experiment is conducted on the breadboard. To begin with, system (1) can be normalized after differential and integral transformation, as follows, where τ 0 is a time constant determined by R 0 and C 0 in the circuit. Considering the timescale transformation, the expressions of the actual excitation signal source are w 1 = Acos(2π f 1 τ 0 t) and w 2 = Acos(2π f 1 τ 0 t), respectively. In other words, the excitation frequency of analog circuit experiment is τ 0 times the result obtained by the numerical analysis. To increase the frequency of the excitation in analog circuit, τ 0 is set to 10000.
Based on (23), an analog circuit is designed to realize the circuit system by using integrated circuits and dis-crete components. The circuit diagram and the experiment equipment are shown in Fig. 14. Multiple UA741, AD633 are used to realize integration, multiplication and addition of Vx, Vy and Vz. Besides, two cosine signals genarated by signal genarator (YB1615P).
In order to obtain the specific component parameters, based on Fig. 14a, the circuit expression also should be deduced as follows Comparing (23) and (24), and setting the circuit parameters as those in Table 3. Finally, the results of the simulation are shown in Fig. 15, and the corresponding hardware experimental results are observed by an oscilloscope (Tektronix MSO3032) and shown in Fig. 16. Consequantly, all of experiment results are consistent with the numerical analysis and simulation results.

Conclusion
In this work, we concentrated on studying a class of multi-timescale non-autonomous dynamical systems with parametric and external excitations. Theoretical analysis shows that the stability of equilibrium will be changed when the parametric and external excitations are introduced. As a consequence, the phenomenon of multi-bifurcation cascade will be formed in the system. Moreover, the numerical analysis results indicated that various kinds of novel multi-bifurcation cascaded periodic, quasi-periodic and chaotic bursting oscillations are observed by changing the frequency of two excitations and the related system parameters. Besides, their formation mechanism is analyzed in detail. Finally, we design an analog circuit by utilizing discrete components to realize the non-autonomous system, and these novel dynamics are verified physically. As a kind of common neurodynamics, the study of the mechanism of bursting oscillation is helpful for the understanding the firing principle of neurons and then to construct spike neural networks that can produce similar dynamics. The study of complex bursting oscillations in this paper is beneficial for some practical engineering applications, such as the establishment of spike neural network.