Multistability and Formation of Spiral Waves in aFractional-OrderMemristor-BasedHyperchaoticLüSystemwith No Equilibrium Points

Multistablity analysis and formation of spiral wave in the fractional-order nonlinear systems is a recent hot topic. In this paper, dynamics, coexisting attractors, complexity, and synchronization of the fractional-order memristor-based hyperchaotic Lü system are investigated numerically by means of bifurcation diagram, Lyapunov exponents (LEs), chaos diagram, and sample entropy (SampEn) algorithm.+e results show that the system has rich dynamics and high complexity. Meanwhile, coexisting attractors in the system are observed and hidden dynamics are illustrated by changing the initial conditions. Finally, the network based on the system is built, and the emergence of spiral waves is investigated and chimera states are observed.


Introduction
A new two-terminal circuit element called memristor characterized by a relationship between the charge and the flux was postulated as the fourth basic circuit element by Chua [1] in 1971. Memristor is short for memory resistor, and its memristance depends on how much electric charge flows through it in a particular direction [2]. Until 2008, researchers at HP Labs firstly published a paper in Nature and reported the realization of the memristor [3]. Since then, applications of the memristor have attracted lots of attention of scholars and a big progress in the fields such as secure communications [4,5], neural network [6,7], chaos [8,9], and simulation of neural network behaviour [10,11] and the memristor oscillators' behaviour [12,13] has been made.
Memristor is a nonlinear electronic element. When the memristor is introduced, the chaotic oscillation of the circuit could be more complex due to its nonlinearity. So the research of memristor-based chaotic systems and its dynamic characteristics has become a new hot spot. In recent years, relevant research results have been obtained [9,[14][15][16]. Fractional calculus is a topic which is more than 300 years old [17]. Lots of physical phenomena show that the fractional-order model can describe the characteristics of the actual system more accurately than the classical integerorder model [18]. Based on this consideration, the application of fractional-order system has attracted more and more researchers' attention [19,20]. When the fractional calculus is introduced into the memristor-based chaotic circuit, a more accurate model for the real situation is obtained. Nowadays, many research studies on chaotic dynamics of the fractional-order memristor-based chaotic system have been conducted [21][22][23][24][25]. Among these studies, Rajagopal et al. [24] analyzed the dynamic properties of the fractional-order memristor system with no equilibrium points through Lyapunov exponents and phase diagrams for the first time. Prakash et al. [25] found rich chaotic behaviours of a novel fractional-order memristor-based chaotic jerk system with no equilibrium points by bifurcation diagram, attractors, and instantaneous phase plots. As a result, the analysis of chaotic dynamical behaviours of the fractional-order memristor-based chaotic systems without equilibrium points is an interesting topic.
Recently, Qiao et al. [26] designed a novel improved memristor-based hyperchaotic Lü circuit without any equilibrium points by adding a linear feedback term, a constant term, and a memristor to the classical Lü chaotic system [27]. e coexisting hidden dynamical behaviours of this system are analyzed. However, the fractional-order chaotic dynamics of the system have not been investigated. Moreover, collective behaviours in the coupled nonlinear network have aroused much research interests, and these networks usually consist of nonlinear chaotic systems [28,29], neural models [30,31], and neural models with memristors [32,33]. In fact, these networks can present different kinds of dynamical behaviours such as synchronization and chimera states. However, there are few reports on the formation of spiral waves in the network of fractionalorder memristor-based hyperchaotic systems with no equilibrium points. us, the dynamics in the network of the new fractional-order memristor-based Lü system need further investigations. e rest of this paper is organized as follows. Section 2 presents the details of the solution algorithm and mathematical model of the fractional-order memristor-based hyperchaotic Lü system. Section 3 investigates dynamics properties and complexity of the system by using bifurcation, Lyapunov exponents, and SampEn algorithm. e multiple coexisting hidden attractors are also illustrated with different initial conditions. Section 4 investigates the synchronization of the system by the ring network structure of the neuron network. Section 5 gives the conclusion.

e Mathematical
Model. e most commonly used memristor-based dynamical system [16] is described as where u and i are the input voltage variable and the input current variable of the memristor, respectively, φ is the internal magnetic flux, α and β are two positive constants, and α + βφ 2 is a memductance function describing the fluxdependent rate of change of charge. Qiao et al. [26] designed the memristor-based hyperchaotic Lü system, and its equations are given by Here, by introducing the fractional calculus to this system, a novel fractional-order memristor-based hyperchaotic Lü system is constructed, and the mathematical model is given by where D q t 0 is the q-order Caputo differential operator [34], x, y, z, and w are state variables, a, b, c, d, and e are the control parameters of classical Lü system [27], and g is the gain of the memristor.
e Caputo fractional-order definition is given in Definition 1.
e Caputo fractional-order derivative definition is given by where q ∈ R + and Γ(·) is the Gamma function.
In order to analyze the dynamic and complexity characteristics of this fractional-order memristor-based hyperchaotic Lü system, we fix the system parameters as a � 36, b � 20, c � 3, d � 5, e � 0.1, α � 4, and β � 0.18 and select the gain parameter g as the only adjustable parameter of the system. It can be obtained by (3): where ∇ q V < 0 means that this fractional-order memristorbased hyperchaotic Lü system is dissipative. It also indicates that the improvement of Lü system and the introduction of memristor and fractional calculus do not change the dissipation proprieties of the original classical Lü system.
has no real solutions and the fractional-order memristor-based hyperchaotic Lü system has no equilibrium. As shown in Ref [26], the integer-order hyperchaotic Lü system has hidden attractors. According to the definition of the hidden attractors [35], the generated attractors and limit cycles of system (3) are all hidden because this fractional-order system investigated in this paper also has no equilibrium points.

e Numerical Analysis Methods.
e predictor-corrector method is an effective method for solving fractionalorder equations. According to Kai et al. [36], it is also called the Adams-Bashforth-Moulton algorithm. Suppose that the fractional-order nonlinear system is defined as 0 is the initial condition of the system, 0 < q ≤ 1, ⌈ · ⌉ is the ceil function, and D q t 0 x(t) is the Caputo derivative. Formula (6) is equivalent to the Volterra integral equation [33,37]: Let h � T/N, t j � jh (j � 0, 1, . . ., N ∈ Z + ), and T be the upper limit for solving this integral equation. en, the correction formula of formula (7) is in which e utilized source Matlab code of the predictor-corrector method can be downloaded online (http://www. sourcecodeonline.com/details/predictorcorrector_pece_ method_for_fractional_differential_equations.html); the function name is "FDE12.m." In [38], the method about how to use "FDE12.m" to solve fraction-order system is described in detail. e LEs (Lyapunov exponents) [39] of system (6) can be calculated as follows. Firstly, suppose Λ k (t) is the eigenvalue of f(x, x(t)), k � 1, 2, . . . , d, where d is the dimension of the system. e orbit x(t) can be given in terms of the tangent linear extension: where v is a vector evolving in the tangent bundle of the manifold and the superscript "T" denotes the transpose. en, the LEs can be defined as the limit eigenvalues of the symmetric operator: where A (t) is the exponential of f(x, x(t)) computed along the trajectory x(t). us, we can obtain e LEs of the trajectory x(t) obtained by system (6) can be determined by

Bifurcations and Lyapunov
Exponents. Dynamics of the system with the variation of g and q are analyzed by means of bifurcation diagram, LEs, and phase diagrams. ree cases are presented.
Case 1. Fix q � 0.998: let the gain parameter g vary from 0 to 20 with step size of 0.04 and the initial condition be [x 0 ; y 0 ; z 0 ; w 0 ] � [1; 0; 1; 0]. Bifurcation diagram and LEs are illustrated in Figure 1. It is shown in Figure 1 that  Figure 2. It is shown in Figure 2 that the system is hyperchaotic when q > 0.972. Meanwhile, when q ∈ [0.965, 0.972], the system vibrates between chaotic and hyperchaotic states. For the rest values of q, the system is convergent. Case 3: Vary parameter g from 1 to 16 with step size of 0.04 and increase derivative order q from 0.94 to 1 with step size 1.2 × 10 − 4 . us, the parameter space is divided as a 400 × 500 grid. As with above cases, the initial condition is given by [x 0 ; y 0 ; z 0 ; w 0 ] � [1; 0; 1; 0]. Maximum Lyapunov exponent-based contour plot in the parameter plane g-q is shown in Figure 3. Chaotic region is observed in a triangular area located in the top left of the q-g parameter plane. When the derivative order q takes different values, the system has relative larger maximum Lyapunov exponents with smaller g. Meanwhile, the maximum Lyapunov exponents of the system decrease with the increase of derivative order q. Finally, it can be found out that the minimum order for chaos is q � 0.927.
Phase portraits of the system with the initial condition [x 0 ; y 0 ; z 0 ; w 0 ] � [1; 0; 1; 0] and different g and q are plotted in Figure 4. Four sets of parameters are given, and they are q � 0.995 and g � 0.3 for Figure 4(a), q � 0.995 and g � 1 for Figure 4(b), q � 0.995 and g � 3 for Figure 4(c), and q � 0.995 and g � 14 for Figure 4(d). Different states including periodic circle, chaotic attractor, and hyperchaotic attractor are observed in the system. According to the above analyses, the conclusion can be drawn that the fractional-order memristorbased hyperchaotic Lü system has rich dynamical behaviours.

SampEn Complexity Analysis.
Sample entropy (SampEn) [40,41] is a measure of the probability size of a new model generated by a time series when its dimension changes. e greater the complexity of the time series is, the greater the probability for a new model that the time series will produce and the greater the entropy holds. e SampEn algorithm can be briefly described as follows: Step 1. Given a time series {x (i), i � 0, 1, 2, . . ., N − 1}, if m represents embedding dimension, then the m-dimensional vectors can be expressed as Step 2. Define the Euclidean distance between X (i) and X (j) as the maximum absolute difference between their corresponding scalar elements, i.e., Step 3. Define the criterion of similarity r as a positive real number. Count the number of d[X(i), X(j)] < r and its ratio to total distance N − m denoted as B r m (i), i � 0, 1, 2, . . ., N − m, j � 0, 1, 2, . . ., N − m, and i ≠ j.
en, the average value of B m (i) can be calculated as Step 4. Similarly, change m to m + 1 and repeat Steps 1 to 3; B (m+1) (r) also can be obtained as Finally, theoretically, the SampEn is defined as e value of SampEn is related to the values of m and r. Usually, m ranges from 2 to 5, and r takes 0.1 to 0.25 times of the standard deviation (SD) of the original time series. is article takes m � 2 and r � 0.2SD.
SampEn complexity of the fractional-order memristorbased hyperchaotic Lü system with parameter g and q is analyzed, and the results are shown in Figure 5. Time series z is used for complexity measure and length of the time series is 10000 with the first 3000 points of data removed. e parameters used for Figure 5(a) are the same as that used for Figure 1. And the parameters used for Figure 5(b) are the same as that used for Figure 2. From Figure 5, we can see that higher complexity is found when the system is chaotic or  hyperchaotic. is result agrees well with the analysis results by the corresponding maximum Lyapunov exponents. Comparing Figure 5(a) with Figure 1(b), it can be seen that the difference of the SampEn complexity is more obvious in Figure 5(a) than the difference of corresponding maximum Lyapunov exponents in Figure 1(b), which indicates that SampEn complexity is more convenient in distinguishing chaotic with nonchaotic states of the system. Let parameter g vary from 0 to 16 with step size of 0.004, derivative order q vary from 0.94 to 1 with step size 1.2 × 10 − 4 , and the initial condition be [x 0 ; y 0 ; z 0 ; w 0 ] � [1; 0; 1; 0]. SampEn complexity in the parameter plane g-q is calculated and shown in Figure 6. Obviously, the higher complexity region in Figure 6 matches well with the larger maximum Lyapunov exponents area in Figure 3. However, unlike the maximum Lyapunov exponent plane, the SampEn complexity plane does not change gradually but is divided into three regions: chaos, periodic cycle, and convergence. e upper right corner in Figure 6 is periodic region. So, the SampEn complexity analysis results are more intuitive and obvious.

Coexisting Hidden Dynamics.
Coexisting attractors can be observed in this system. Fix q � 0.998 and g � 5; four phase diagrams are shown in different forms in Figure 7 because of the different initial conditions. When x 0 � y 0 � z 0 � w 0 � 0 and x 0 � y 0 � z 0 � w 0 � 1, the system is periodic but has different number of cycles. When x 0 � y 0 � z 0 � w 0 � 3 and x 0 � y 0 � z 0 � w 0 � 5, the system is chaotic, but Figure 7(d) shows the weak chaotic state. When g varies from 4.8 to 5, the bifurcation diagrams for the corresponding 4 different initial conditions are shown in Figure 8. e blue points are for x 0 � y 0 � z 0 � w 0 � 0, red points are for x 0 � y 0 � z 0 � w 0 � 1, blue points are for x 0 � y 0 � z 0 � w 0 � 3, and magenta points are for x 0 � y 0 � z 0 � w 0 � 5. It can be seem from Figure 8 that the system trajectory shows the multiple stable phenomena of coexistence for different periodic limit cycles, chaos, and weak chaos due to the different initial conditions.
In order to further study the effect of initial conditions on the dynamics of the system, the SampEn complexity chaos diagram for different parameters is illustrated in Figure 9. Comparing Figure 9(a) with Figure 9(b), it can be seen that the effect of initial conditions on system dynamics is multidimensional. Each initial component will play a role in the dynamics of the system. Comparing Figure 9(c) with Figure 9(d), we can find that when the derivative order q is different, the effect of initial value on system dynamics is also different. By comparing the four figures in Figure 9, it can be found that the SampEn complexity diagrams are different with different initial conditions and SampEn algorithm is an effective method for the analysis of the hidden dynamics in the system.
Compared with the basin attraction plots, the SampEnbased chaos diagrams illustrate the dynamics and complexity of the system with different initial conditions. If the complexity is different, the state of the system is different.
us, it can show the multistability of the system in the initial plane through a complexity of view. Since complexity Mathematical Problems in Engineering

e Structure of the Ring Network.
To further capture the complex dynamics of the system, we construct a ring network of nonlinear systems which are coupled locally to the nearby 3 nodes with a coupling strength of D. e ring network structure of the neuron network is depicted in Figure 10 and can be numerically defined as where i = 1 to N. N is the number of nodes. e set nodes contain the three nodes near the target node. As shown in Figure 10, the nearby nodes of node 1 are node 2, node 3, and node 49, while the nearby nodes of node 36 are node 34, node 35, and node 38. To find the nodes, Algorithm 1 is used. In order to derive the dynamics of the neuronal ring network, we also fix the system parameters as a � 36, b � 20, c � 3, d � 5, e � 0.1, α � 4, and β � 0.18 and vary the gain parameter g. In Figure 10, the built network has 50 nodes, and thus N � 50.
; end if ALGORITHM 1: e algorithm to find the nodes near node i.

Numerical Simulation of the Fractional-Order Network.
Spatiotemporal patterns of the model with the variation of g and q are analyzed for different coupling strengths. ree cases are presented. Case 1. Fix q = 0.998, g � 3, and various coupling strengths.
e spatiotemporal patterns are shown in Figure 11. In these figures, the level of excitability of the network is illustrated by colours. We can see that when the coupling strength D = 2 and 4, only unsynchronized nodes are shown in Figures 11(a) and 11(b). When the coupling strength D increases to 6, some synchronized neurons are found in Figure 11(c), but there are still    Figure 13 shows that the systems are not synchronized for the coupling strength D = 1. With the coupling strength increasing gradually, the number of synchronized systems increases. When the coupling strength D = 6, all nodes are in a local synchronization state between different nodes.
As shown in Figures 11-13, the behaviours of the network are different with different coupling strengths. When the coupling strength is small, the network does not show synchronization. For instance, it is shown in Figures 11(a) and 11(b) that the network is not synchronized. When the coupling strength is large, the network can be synchronized. For example, Figure 12(d) shows that there exists synchronization. e chimera state means coexistence of synchronization and asynchronization in the network. According to Figures 11-13, there are different kinds of chimera state in the network. e network with different fractional-order derivative order has different dynamics. As for the coupling strength for chimera states in the network, it should be proper. According to the analysis results, when the coupling strength D ∈ [2,4], the chimera states are observed. Moreover, the derivative order also affects the dynamics of the network. In Case 1, we fix q � 0.998 and in Case 3 we let q � 0.97, and it is shown in Figures 11 and 13 that the dynamics of the network are different. It is easier for the network with q � 0.97 to synchronize or to show chimera states.

Conclusions
In this paper, a fractional-order memristor-based hyperchaotic system with no equilibrium points has been numerically analyzed by employing the predictor-corrector method. It shows that the system has rich dynamical behaviours, and different states including periodic circle, chaotic attractor, and hyperchaotic attractor are plotted in the system. e minimum derivative order of generating chaos is found, and it is q � 0.927. Meanwhile, through numerical simulation, the coexisting hidden attractors are observed with different initial conditions. Moreover, the SampEn complexity verifies the hidden dynamic characteristics of the system and it shows that this complexity analysis algorithm can be employed to detect the multiple coexisting attractors. Finally, by constructing the ring network structure of the system, it can be found that with increase of the coupling strength, the state of the network is different. e network can have the global synchronization state with large coupling strength, while in most cases, the network shows chimera states with proper coupling strength.

Data Availability
All data, models, and code generated or used during the study appear in the submitted article.

10
Mathematical Problems in Engineering 2019M652791), and the Postdoctoral Creative Talent Support Programme (grant no. BX20180386).