Dynamics of stimuli-based fractional-order memristor-coupled tabu learning two-neuron model and its engineering applications

External stimulus has an impact on the functional behavior of the biological nervous system, and appropriate stimulus helps the organism to maintain neural function. Inspired by this, the effects of different external stimuli on the dynamical behaviors of neuron model are studied in this paper. Firstly, a fractional-order (FO) memristor-coupled tabu learning two-neuron model whose equilibrium points are symmetric about the origin and unstable is constructed. The dynamical behaviors of this neuron model under different stimuli are further discussed, which are no external stimulus, external forced current stimulus, and electromagnetic radiation (EMR), respectively. The neuron model without external stimulus has periodic attractors and transient chaos, but applying external stimulus to one of the neurons can produce abundant two-scroll chaotic attractors and multistability; especially, when the neuron model is stimulated by EMR, it can generate hyperchaotic attractors that have not been observed in the tabu learning neuron model before. Besides, the transient transition behaviors of the model under different stimuli are also studied. Then, a pseudo-random number generator is designed and its random performance is tested with NIST suite. Finally, it is applied to voice encryption, and the result shows that it has good encryption effect. Therefore, it can be said that the FO memristor-coupled tabu learning two-neuron model has superior randomness and is suitable for chaotic-based engineering applications.


Introduction
Biological nervous system has been studied in order to explore the mechanism of artificial intelligence in recent years [1,2]. The nervous system is composed of a large number of interconnected neurons, which exchange information by receiving, integrating, transmitting and exporting energy [3,4]. In order to simulate the dynamical behavior of neurons, many artificial neurons and neural network models have been developed and utilized. Among them, tabu learning neural network proposed by Beyer and Ogier is a neural network model that uses tabu search algorithm to ensure diversified and effective exploration to achieve global optimization [5,6]. There are a few studies on the complex dynamical behavior of tabu learning neuron model. It is necessary to further study its complex discharge dynamics, which has important scientific significance and application value.
A large number of physical and biological experiments have confirmed that many neurological diseases such as Alzheimer and epilepsy [7] are caused by abnormal discharge of neurons [8,9], and appropriate stimulus helps to reduce the symptoms of psychiatric patients, which is conducive to their health [10]. Hence, studying the firing behavior of neurons under external stimulus plays a positive role in understanding the pathogenesis of neuronal diseases and then contributes to the treatment and prevention of such neurological diseases. It has been proved that the dynamics of a neuron or neural network can be affected by external stimulus, such as noise [11], electric field [12], and magnetic field. At the same time, EMR has become a factor that cannot be ignored because it is almost full of human living environment. By introducing an additional membrane current into the Hodgkin-Huxley neuron model, Li et al. [13] first found that EMR can inhibit the firing activity of a single neuron and regulate the collective electric activities of neuronal network. In order to explore the influence of EMR distribution on the nonlinear dynamics of a neural network with n neurons, Lin et al. [14] investigated the chaotic dynamics of neural network model consisted of three neurons by stimulating different number of neurons. Ma et al. [15] confirmed that electric field stimulation can cause distinct mode transition in electrical activities of neuron exposed to different kinds of electric field. Ge et al. [16] constructed the improved HR neurons in the feed-forward multilayer networks and studied the effects of EMR, synaptic weight, and noise intensity on the propagation of the subthreshold excitatory postsynaptic current signal and the neural synchronization. The dynamics of tabu learning neuron model stimulated by different disturbances, namely external forced current stimulus, EMR, and both EMR and external forced current stimuli, are discussed by Li et al. [17].
Memristor applied in neural network can not only be used as synapse or autapse to form memristive synaptic neural network model [18], but also characterize EMR by describing the relationship between magnetic flux and membrane potential to form a neural network model under EMR [19]. What is more surprising is that the memristive neural network model has more complex dynamical characteristics than the ordinary neural network model, which has been proved in many studies [20][21][22][23]. Taking the memristor as a self-feedback synapse of a neuron, Lai et al. [24] present a design of a new Hopfield neural network that can yield multi-double-scroll attractors and can effectively control the number of double scrolls contained in an attractor. Hou et al. [25] proposed a memristive self-connection synaptic weight-based tabu learning neuron model to study the coexisting infinitely many nonchaotic attractors composed of mono-periodic, multiperiodic, and quasi-periodic orbits. By using memristor to characterize EMR, Ren et al. [26] proved that EMR can not only contribute to the discharge of neurons under positive feedback coupling, but also enhance synchronization behaviors of coupled neurons under negative feedback coupling. The initialinduced extreme multistability of heterogeneous memristive bi-neuron network and the synchronous firing activities in homogeneous memristive bi-neuron network are investigated by Bao et al. [27] through bidirectionally coupling two three-dimensional heterogeneous or homogeneous Morris-Lecar neurons with a memristor synapse. Lin et al. [28] established a memristive Hopfield neural network by utilizing the multi-stable memristor as a memristor synapse, which generated infinitely many coexisting chaotic attractors that enjoy the same structure.
Multistability and extreme multistability are important topics in the study of memristive nervous system. Multistability refers to a variety of steady-state behaviors caused by the change of initial state when the system parameters are unchanged, that is, a variety of attractors exist in the system under different initial conditions [29]. In particular, when an infinite number of attractors exist at the same time, this phenomenon is defined as extreme multistability [30]. Therefore, extreme multistability is a special case of multistability. If the coexisting multiple attractors have different positions but the same shape, it is homogeneous multistability; otherwise, it can be considered as heterogeneous multistability [31]. The memristive nervous system with multistability characteristics has great flexibility to switch from one stable state to another, which better reflects the brain's firing behavior, so it has attracted extensive attention from researchers. By introducing the nonvolatile locally active memristor into the neural network model composed of three Hopfield neurons, Li et al. [32] revealed the coexistence phenomenon of multiple stable modes. Considering the influence of electromagnetic induction and external stimulus, extreme multistability in the continuous non-autonomous memristive Rulkov model was discovered by Xu et al. [33]. A new no-equilibrium HR neuron model with memristive electromagnetic induction is proposed by Zhang et al. [34], which breeds the interesting phenomenon of hidden homogeneous extreme multistability. Bao et al. [35] present an improved nonautonomous memristive FitzHugh-Nagumo circuit, from which coexisting infinitely many attractors are obtained.
Compared with integer-order calculus, FO calculus has the advantages of infinite storage and greater degrees of freedom [36] and can more accurately describe the mathematical models in bioengineering, signal processing, chaos theory, and other fields [23,37]. In recent years, FO calculus has become a research hotspot of memristive neural network. On the one hand, memristor has nonvolatile and memory characteristics, which matches the genetic and memory characteristics of FO calculus in the description process [38]. On the other hand, the action potential or membrane potential of neurons is caused by the directional movement of ions, and the attenuation rate of ion current is related to the FO derivative [39]. Therefore, it is reasonable and accurate to use FO calculus to simulate memristive neural networks. Generally speaking, the firing behaviors of FO neuron model are more complex than that of integer-order neuron model. Yu et al. [40] remodeled the slow ion channel of FO memristive HR model with a linear FO memristor, took the FO memristive HR model as the small perturbation system of fast subsystem, and found that the system showed different membrane potentials for different orders. Based on the FO memristive synaptic-coupled neuron model, Xin et al. [41] discussed its synchronization behaviors and found that the synchronization transition was affected by the fractional order. By using the coupled FO locally active memristor to simulate the synaptic crosstalk in Hopfield neural network, Ding et al. [42] found that the HNN model with coupled locally active memristor has multi-stability under different fractional order and coupling coefficient. Among many memristive neuron models, there is no research on memristive tabu learning neuron model in the field of FO calculus. Therefore, it is necessary to study the dynamical behaviors of FO memristive tabu learning neuron model. Inspired by the above discussion, the dynamical research and engineering application of memristive tabu learning neuron model are further expanded. Firstly, a new FO memristor-coupled tabu learning two-neuron model is proposed in this paper. The chaotic dynamics of the neuron model under different stimuli are further studied, namely no external stimulus, external forced current stimulus, and EMR. Its main contributions are summarized as follows: (1) There are a few studies on introducing the concept of memristor simulating synapse into tabu learning neurons. Further, the study of unidirectional coupling of two tabu learning neuron models through memristor synapse has not been reported. Therefore, based on the properties of fractional order and memristor, we coupled two adjacent tabu learning neurons unidirectionally through memristor synapse, established a FO memristor-coupled tabu learning two-neuron model, and further theoretically analyzed its equilibrium point and stability. (2) Based on the FO memristor-coupled tabu learning two-neuron model, the dynamical behaviors of the neuron model under different stimuli are further studied, that is, no external stimulus, external forced current stimulus, and EMR. It is found that the neuron model without external stimulus has periodic attractors and transient chaos, but abundant two-scroll chaotic attractors can be generated by applying external stimulus to the first neuron. Therefore, the system shows more complex firing behaviors in the presence of external stimuli. (3) We also note that the firing behaviors of the neuron model under the stimulation of EMR are much more complex than that under the stimulation of external forced current. Under the stimulation of EMR, the neuron model has multistability, hyperchaotic firing behavior, and multiple transient transition behaviors. This is of great research significance for the treatment and prevention of neurological diseases by applying appropriate stimulation. (4) The FO memristor-coupled tabu learning twoneuron model is applied to the PRNG. The randomness of the two types of attractors is detected and compared by NIST suite. It is further applied to voice encryption, and the statistical results of encryption are given.
The rest of this paper is arranged as follows. In Sect. 2, a FO memristor-coupled tabu learning twoneuron model is proposed, and the stability of the equilibrium point is analyzed. In Sect. 3, the complex dynamical behaviors of neuron model under three different external stimuli are studied by using bifurcation diagrams, Lyapunov exponents (LEs), and local attraction basins. In Sect. 4, the neuron model is applied to the PRNG. Section 5 introduces the application in voice encryption. The last section summarizes the whole paper.

Tabu learning two-neuron model
Hopfield neural network is usually used to minimize an objective function. When dealing with optimization problems, Hopfield neural network model can be described as [43,44]: where V i ¼ f u i ð Þ, f Á ð Þ is the activation function, u i is the state variable of the i-th neuron, C i and R i are positive constants representing capacitance and resistance, T ij is the connection weight from the j-th neuron to i-th neuron, I i represents the input current to the i-th neuron.
When the linear similarity function is selected, the equation of state of the tabu learning is: J i satisfies the following: where c and d are the positive constants that represent the memory decay rate and the learning rate. Combining Eqs. (2) and (3), a tabu learning twoneuron model can be written as: T is the state variable of the neurons, y i ¼ y 1 ; y 2 ½ T is the tabu learning state variable, b ij represents a synaptic weight matrix 2 Â 2, which describes the connection strength between the two neurons, a i ¼ 1=R i , c i and d i are positive constants, T is the activation function, and T is the externally input current.

Fractional calculus
There are three commonly used definitions of FO derivatives, namely the Grünwald-Letnikov (GL), Riemann-Liouville (RL), and Caputo definitions. In mathematical theory, the RL derivative is more commonly used than the Caputo derivative, but in physical systems, the Caputo derivative is more practical than the R-L derivative. Since the Caputo definition allows the use of initial values of classical integer-order derivatives with clear physical interpretation, Caputo derivative is usually utilized in engineering applications. Therefore, the Caputo derivative will be used to model the system in this paper. For the function f t ð Þ, the q-order derivative under the Caputo definition is defined as [45,46]: where D q t 0 is the fractional derivative of order q 2 R þ with respect to time t ð Þ, Cðm À qÞ is the gamma function of m À q formulated as Cðm À qÞ ¼ R þ1 0 e Àt t mÀqÀ1 dt, t 0 is the initial time and m is an integer satisfying m À 1\q\m.

Model description
Synaptic coupling plays a positive role in signal transmission between neurons. When two adjacent neurons are coupled through synapses, the electromagnetic induction current will be generated if there is a membrane potential difference between the two neurons, which can be characterized by a fluxcontrolled memristor synapse. The memristor synapse is expressed as [47]: where v m ¼ x À u, i m , u represent the membrane potential difference, induced current, and inner flux state variable triggered by the membrane potential difference of two coupling neurons, respectively. a, b, g, l are the four internal positive parameters of the memristor, which are set to a ¼ 1: When the order is selected as q ¼ 0:9, the driving signal is a sinusoidal voltage source where v m and f are the amplitude and frequency. The pinched hysteresis loops of the memristor are given.
When v m ¼ 1:7 V is maintained unchanged and f is selected as 500 Hz, 1 KHz and 2 KHz, respectively, the v À i curve is shown in Fig. 1a, while f ¼ 1 KHz is fixed, and v m is assigned as 1.6 V, 1.7 V, 1.8 V, the v À i curve is shown in Fig. 1b. It can be seen that the v À i curve is the hysteresis loop pinched at the origin, and its side lobe area decreases with the increase in frequency and finally shrinks as a linear function. Therefore, the memristor (7) meets the basic characteristics of memristor. Tabu learning neuron model (4) requires the activation function to be bounded and differentiable. Therefore, the nonlinear activation function is selected Based on the model of memristor synapse in (7), a FO memristor-coupled tabu learning two-neuron model is derived, which is expressed as: where q 2 0; 1 ð , x and u stands for the neuron state variables. y and v correspond to the tabu learning state variables; the parameter r represents the coupling strength of electromagnetic induction between the two neurons. Furthermore, neuron internal constants a 1 ; a 2 ; c 1 ; c 2 ; d 1 ; d 2 are selected as a 1 ¼ 2; The simplified topology of FO memristor-coupled tabu learning twoneuron model is shown in Fig. 2.

Equilibrium point and stability analysis
The stability of equilibrium point is an important property of nonlinear dynamical system. The equilibrium point of FO memristor-coupled tabu learning two-neuron model, which is denoted by S ¼ x; y; u; v; u ð Þ , can be solved by the following equations: x; u ð Þ can be seen as the intersection of the curves drawn by Eqs. (9a) and (9b); therefore, the values of x and u can be solved by the graphical analysis method. When q ¼ 0:9, the coupling strength r ¼ 0:7, other parameters remain unchanged, and Eq. (9a) and (9b) is plotted in Fig. 3. The solution of x and u is the intersection of two functions in Fig. 3. Substituting x and u into Eq. (10), we get five equilibrium points,p 0 ¼ 0; 0; 0; 0; 0 ð Þ ; p 1 ¼ 0:2854; À0:8337; ð À0:7592; 1:9217; À1:9182Þ;   The Jacobian matrix of the system at the equilibrium point S ¼ x; y; u; v; u ð Þis shown below: where n ¼ a À bu 2 .
The eigenvalues must satisfy the following formula for the equilibrium point of FO model to be asymptotically stable.
In other words, when q satisfies Eq. (14), the equilibrium point is unstable.
For the point p 0 ¼ 0; 0; 0; 0; 0 ð Þ , the corresponding eigenvalues can be calculated as À1:2000, 4.7733, À5:2842, À0:2155, À0:8135, which indicates the equilibrium point p 0 is unstable. Similarly, the other four equilibrium points in this paper are unstable. As a result, the unstable equilibrium point may lead to the occurrence of chaos or periodic oscillations in the FO memristor-coupled tabu learning two-neuron model.

Dynamics of the neuron without external stimulus
The formula of the model without external stimulus is the same as Eq. (8). The connection topology is shown in Fig. 4. The dynamical phenomena of the model are researched by using bifurcation diagrams, LEs, phase diagrams, and time-domain waveforms. The bifurcation diagram depicts the characteristics of the system performance varying with the system parameters; the value of the LEs represents the average divergence or convergence speed of adjacent orbits in the attractor along a certain direction. Set And set the initial values as 1; 0; 1; À1; 1 ð Þ , À1; 0; À1; 1; À1 ð Þ , respectively. The coupling strength r changes within the range of 0:2; 1:5 ½ , the bifurcation diagrams, and LEs corresponding to r are  Fig. 5. Figure 5 shows that the system basically presents the same periodic firing pattern under the two groups of initial values. Meanwhile, due to the existence of transient chaos, the maximum Lyapunov exponent in some places is slightly greater than zero, but it is basically consistent with the bifurcation diagrams. The three-dimensional phase diagram and the corresponding time-domain waveform are shown in Fig. 6. The phase diagram reflects the change of the system state and can be used to directly judge the simple dynamical behaviors of the system. In Fig. 6a, c, the red-, blue-period-4 attractors correspond to the initial value AE1; 0; AE 1; Ç 1; AE 1 ð Þwhen r ¼ 0:7. Figure 6b, d shows that when the coupling strength r ¼ 1:4, and the initial value is AE1; 0; AE 1; Ç 1; AE 1 ð Þ ; the system is in the state of red-, blue-period-2. It is worth mentioning that Model (8) produces symmetric periodic attractors under the coordinate transformation x; y; u; v; u ð Þ$ À x; Ày; Àu; Àv; Àu ð Þ . The emergence of transient chaos is the result of system evolution with time, which is reflected in the transition from transient chaos to long-term stable periodic state. The time-domain waveform and threedimensional phase diagram of Fig. 7 explore the transient chaotic behavior of Model (8) without external stimulus. In order to better explain the transient behavior of the system, the time-domain diagram of variable u is given, and the two-dimensional phase diagram is supplemented.
Coupling strength r is selected as r ¼ 0:7. On the one hand, the initial value is set to 1; 0; 1; À1; 1 ð Þ ; when t 2 0; 94:2 ð Þ, the system is in a chaotic state and is marked as green in the time-domain waveform and phase diagram in Fig. 7. The system is in the period-4 state when t 2 94:2; 500 ð Þ , which is marked red. On the other hand, the initial value is À1; 0; À1; 1; À1 ð Þ , when t 2 0; 94:2 ð Þ, the system is in magenta chaotic firing pattern, and the system is in blue-period-4 firing pattern when t 2 94:2; 500 ð Þ .

Dynamics of the neuron under external forcing current
This section considers that the tabu learning neuron 1 is influenced by external forcing current, while tabu learning neuron 2 has no external stimulus. The effects of the frequency F of sine pulse and the coupling strength r on the dynamical behaviors of the FO memristor-coupled tabu learning two-neuron model are studied. The model of FO memristor-coupled tabu learning two-neuron model under external forcing current can be established as follows: where the external forcing current I 1 ¼ A sin 2pFt ð Þ, I 2 ¼ 0. The connection topology of the Model (15) is shown in Fig. 8.
When q ¼ 0:9, r ¼ 0:7, A ¼ 0:8, and other parameters remain unchanged, the coexisting bifurcation behaviors with the frequency F changes in the range of 1; 2:5 ½ are shown in Fig. 9. The red and blue trajectories represent the initial values of 1; 0; 1; À1; 1 ð Þand À1; 0; À1; 1; À1 ð Þ , respectively. It can be seen that the bifurcation diagrams are basically consistent with the dynamical behaviors of LEs of Fig. 9b. As shown in Fig. 9a, when F is in the range of ð1; 1:43Þ, the system is in the chaotic firing pattern, then enters the period-3 firing when F ¼ 1:43, where the quasi-period-3 firing is in the range of ð1:52; 1:59Þ, until F ¼ 1:8, the system switches to chaotic firing again. Within the range of ð2:35; 2:41Þ, the system is in periodic firing pattern, finally, it enters into chaotic firing when F ¼ 2:41.
The Model (15) generates rich symmetric coexisting firing behaviors, which is well shown in the phase diagrams in Fig. 10, and the time-domain waveforms are completely corresponding to them. As shown in Fig. 10a, d, when F ¼ 1:5, the system is in the symmetric-period-3 firing pattern. When F ¼ 2, the symmetric two-scroll chaotic firing pattern appears in Fig. 10b, e. When F ¼ 2:4, Fig. 10c, f shows that the system shows symmetric-period-4 firing behavior.
When q ¼ 0:9, A ¼ 0:8, and F ¼ 1:5 remain unchanged, the bifurcation behaviors of r in the range of 0:7; 3:5 ½ are shown in Fig. 11. Due to the symmetric coexistence, the bifurcation routes of the two groups of initial values are consistent. Red, blue trajectories start from the period-3 firing, enter the chaotic firing at r ¼ 0:73, return to period-3 firing again at r ¼ 1:05, and then into a the period-6 firing via period-doubling bifurcation at r ¼ 1:4. After that, red, blue trajectories enter the chaotic firing again at r ¼ 1:55 until the end, during which two small period windows appear, for example, when r ¼ 1:85, the system presents the period-3 firing, and the system presents the period-4 firing when r ¼ 3:02. In addition, the LEs correspond well with the bifurcation diagrams. Figure 12 shows the coexistence of different topologies at different coupling strengths r and Þ , and the blue attractor corresponds to the initial value À1; 0; À1; 1; À1 ð Þ . To further explore the effects of the external current frequency F and the coupling strength r on the firing behaviors of neuron model, two-parameter bifurcation diagram of F and r is presented. Two-parameter bifurcation diagram is measured with two parameters, where different types of dynamical behaviors are marked with different colors. As shown in Fig. 13a, when the initial value is 1; 0; 1; À1; 1 ð Þ , the blue region represents the chaotic state, the cyan region represents the period-1 firing, the green region represents the period-2 firing, the yellow region represents the period-3 firing, the orange region represents the period-4 firing, and the red region represents the multiperiod firing with period greater than 4. In addition, Fig. 13b shows a two-dimensional standard Lyapunov exponent diagram on the parameter plane, corresponding to the first Lyapunov exponent, where the red and yellow regions represent chaotic states (k max [ 0) and periodic states (k max \0), respectively. It can be seen that the system has abundant firing behaviors, in which the blue chaos is the main firing pattern.
Under the stimulus of external current, we find that Model (15) produces two new phenomena when the initial value is 1; 0; 1; À1; 1 ð Þ : intermittent chaos and state jump. Time-domain diagram of Fig. 14a shows that local chaos and local non-chaos oscillation transition appear alternately when F ¼ 1:3. The local chaotic state is shown by the cyan chaotic attractor in Fig. 14b, c, and the local non-chaotic state is shown by the blue periodic attractor. Meanwhile, when F ¼ 1:9, Fig. 14d-f shows that the green two-scroll chaotic state jumps to the red two-scroll chaotic state at t ¼ 425.

Dynamics of the neuron under EMR
Sections 3.1, 3.2 shows that the FO memristor-coupled tabu learning two-neuron model presents relatively simple dynamical behaviors without external stimulus and under external forcing current. Therefore, this subsection considers applying EMR to this neuron model and discusses its rich dynamical behaviors.

FO memristor-coupled tabu learning twoneuron model under EMR
In the electrophysiological environment, the electrical activity of neurons may be affected by EMR. From the perspective of electromagnetic induction theory, magnetic flux can be used to describe EMR effects. The i ¼ W q u 0 ð Þv ¼ 2 tanhðku 0 Þ À tanhðku 0 þ 1:5kÞ ð À tanhðku 0 À 1:5kÞÞv; where v, i represent the input voltage, the induced current caused by the change of magnetic flux u 0 , k is the internal parameter of the memristor, and the order q satisfied 0\q\1. W q ðu 0 Þ denotes the memductance function.
When the order is selected as q ¼ 0:9, the inner parameter of the FO flux-controlled memristor is set to k ¼ 2:5; the driving signal is a sinusoidal voltage where v m and f are the amplitude and frequency. The pinched hysteresis loops of the memristor are given. When v m ¼ 1:7 V is maintained unchanged and f is selected as 500 Hz, 1 KHz and 2 KHz, respectively, the v À i curve is shown in Fig. 15a, while f ¼ 1 KHz is fixed, and v m is assigned as 1.4 V, 1.7 V, 2 V, the v À i curve is shown in Fig. 15b. It can be seen that the v À i curve is the hysteresis loop pinched at the origin, and its side lobe area decreases with the increase in frequency and finally shrinks as a linear function. Therefore, the memristor (16) meets the three basic characteristics of memristor [48,49].
As shown in Fig. 16, we think that the tabu learning neuron 1 is affected by EMR; the tabu learning neuron 2 has no external stimulus. The neuron model can be established as follows: where E 1 ¼ qW q u 0 ð Þx ¼ qð2 tanhðku 0 ÞÀ tanhðku 0 þ 1:5kÞ À tanhðku 0 À 1:5kÞÞx, E 2 ¼ 0. u 0 denotes the magnetic flux across cell membrane of the tabu learning neuron 1. The term qð2 tanhðku 0 Þ À tanhðku 0 þ 1:5kÞÀ tanhðku 0 À 1:5kÞÞx is an induced current caused by the change of magnetic flux and field, where q represents the strength of external EMR.

Coexisting bifurcation behaviors
Under the stimulus of EMR, the Model (17) has a variety of firing behaviors. By numerical simulation methods such as bifurcation diagrams, Les, and local attraction basins, the complex coexisting multiple firing behaviors are revealed when the EMR strength q and the coupling strength r change.
When the order q ¼ 0:9, the coupling strength r ¼ 0:7, three sets of initial values are selected, the max-spike bifurcation diagrams and the first two LEs of the EMR strength q in the region 0:8; 6:5 ½ are drawn in Fig. 17a, b. When the initial value is 1; 0; 1; À1; 1; 1 ð Þ , the red trajectory starts from period-4 firing and becomes quasi-periodic firing at q ¼ 1:25, then enters period-1 firing at q ¼ 3:03, jumps to period-3 firing pattern at q ¼ 4:5. With the continuous increase in q, the system enters the chaotic ½ , the dynamical behaviors represented by LEs and bifurcation diagrams are inconsistent in some specific time intervals. When the initial value is À5; À5; 1; À1; 1; 1 ð Þ , the blue trajectory begins with period-4 firing and enters the quasiperiodic firing at q ¼ 1:79. After that, it experiences chaotic firing in a very small range of 3:87; 4:08 ½ and then converts to the periodic firing. It jumps to the quasi-periodic firing when q ¼ 4:71, gradually evolves into the periodic firing, and finally breaks into the multi-periodic firing after q ¼ 5:62. In the end, the purple trajectory starts with the hyperchaotic firing when the initial value is À1; 0; À1; 1; À1; À1 ð Þ , then degenerates into the chaotic firing at q ¼ 1:34, and after that, enters into large-range quasi-periodic firing at q ¼ 1:5. When q increases to 3.87, it becomes chaotic firing again and finally transitions to the periodic firing at q ¼ 4:08. Further, it can be seen from LEs of Fig. 17b that the Lyapunov exponent agrees with the dynamical behavior of the bifurcation diagram.
The EMR strength is fixed as q ¼ 1, other parameters remain unchanged, and the coupling strength r changes in the range of 0:2; 4 ½ . The corresponding bifurcation behaviors are shown in Fig. 18. The trajectories of the three colors correspond to three different initial values: green, red, and blue trajectories represent the initial values À1; 0; À1; 1; À1; À1 ð Þ , À2:3; À2:2; 1; À1; 1; 1 ð Þ , and 1; 0:5; 1; À1; 1; 1 ð Þ , respectively. It can be observed that the dynamical behaviors of the system with r are different under three different initial values. When the initial value is À1; 0; À1; 1; À1; À1 ð Þ , the system starts from hyperchaotic firing, degenerates into the chaotic firing at r ¼ 1:17, is in the multi-periodic firing pattern in the range of 1:34; 1:45 ½ , and then enters the quasiperiodic firing. It can be seen from LEs at the top of Fig. 18b that LE 1 ¼ 0 and LE 2 ¼ 0. Until r ¼ 2:42, the quasi-periodic firing changes into the period-3 firing.
When the initial value is À2:3; À2:2; 1; À1; 1; 1 ð Þ , it can be found by observing the red trajectory that the system starts from the period-2 firing, enters the period-4 firing via perioddoubling bifurcation at r ¼ 0:37, then returns to the period-2 firing by the reverse period-doubling bifurcation at r ¼ 0:89, changes to the quasi-period-3 firing at r ¼ 2:02, and enters period-3 firing at r ¼ 2:34 until the end. Finally, before r ¼ 2:75, the blue trajectory is consistent with the red trajectory. After that, the blue trajectory enters the quasi-periodic firing and then enters the period-2 firing when r ¼ 3:58. Figure 18 shows that when EMR is added, the dynamical behaviors of the coupling strength r are richer than that of the first two cases, indicating that the model has abundant coexistence multi-stability under external electromagnetic stimulation.

Coexisting attractors and multistability
In order to further discuss the multiple coexisting attractors of neuron model under EMR, the local   attraction basins related to the initial conditions are given, in which different types of dynamical states are marked with different colors. When q ¼ 0:9, the initial condition ðuð0Þ; vð0Þ; uð0Þ; u 0 ð0ÞÞ¼ ð1; À1; 1; 1Þ, the coupling strength r ¼ 0:7, Fig. 19 shows the local attraction basins on the initial plane of x 0 ð Þ À y 0 ð Þ when the EMR strength q is q ¼ 5:28 and q ¼ 6:5, respectively, in which seven different colors represent different types of firing behaviors. Here, BP3, RCH, OP2, YQP, CP2, WP2, BP5, RP5, YP2, CP1, and UB represent brown-period-3, red-chaos, orange-period-2, yellow-quasi-period, cyan-period-2, wathet-period-2, brown-period-5, red-period-5, yellow-period-2, cyan-period-1, and unbounded, respectively. At the same time, the coexisting three-dimensional phase diagrams corresponding to the local attraction basins are given. Six groups of initial values are selected in Fig. 20a, b, which correspond to the six bounded states of the local attraction basin in Fig. 19a, b, respectively.
As can be seen from Fig. 19a, when q ¼ 5:28 is fixed, the period-3 attractor is represented by the brown region, the chaotic attractor is represented by the red region, and the brown-period-3 attractor at the top right of Fig. 19a is dotted in the red-chaotic region completely. The orange region represents the period-2 attractor, the yellow region represents the quasi-period attractor, and the orange-period-2 and the yellow quasi-period region dominate the initial plane. Cyanperiod-2 attractor is largely intertwined with wathetperiod-2 attractor. Tiny blue region represents unbound. When q increases to 6.5, r ¼ 0:7 remains unchanged; Fig. 19b shows that the brown-period-5 attractor, the red-period-5 attractor, and the cyanperiod-1 attractor are nested in the upper part of Fig. 19b. The yellow-period-2 attractor and the wathet-period-2 attractor are interleaved in the lower part of Fig. 19b. It is worth noting that the riddled attraction basin is displayed in a large region, which implies that the model in this case is sensitive to the initial conditions.

Hyperchaotic attractors
According to Sect. 3.3.2, when the initial value is selected as À1; 0; À1; 1; À1; À1 ð Þ , the system can produce hyperchaotic attractors with the change of EMR strength q and coupling strength r. Hyperchaotic firing behavior should have at least two LEs greater than zero. Figure 21c shows that the corresponding first two LEs of the system in the region ½ . This means that the system produces hyperchaotic attractors. The phase diagrams of the two-scroll hyperchaotic attractor are shown in Fig. 23 when r ¼ 0:7, q ¼ 1. The corresponding first two LEs are LE 1 ¼ 0:2125, LE 2 ¼ 0:0206, indicating that it is a hyperchaotic attractor.
In order to further explore the influence of parameters on hyperchaotic firing behavior, when the initial value is À1; 0; À1; 1; À1; À1 ð Þ , the two-dimensional bifurcation diagram is plotted in the region q 2 0:8; 6:5 ½ , r 2 0:2; 1:5 ½ . As shown in Fig. 24a, in which wathet region represents periodic firing, cyan region represents quasi-periodic firing, yellow region represents chaotic firing, and orange region embedded in yellow region represents hyperchaotic firing. In  In addition, Fig. 24b, c shows the two-dimensional standard Lyapunov exponent diagrams on the q; r ð Þ parameter plane. Figure 24b corresponds to the first Lyapunov exponent, where red and yellow regions, respectively, indicate that the value of Lyapunov exponent is positive and close to zero. Figure 24c corresponds to the second Lyapunov exponent, where the red, yellow, and blue regions represent the values of the Lyapunov exponent as positive, close to zero and negative, respectively. Hyperchaotic firing occurs when red regions in Fig. 24b, c appear at the same location in the diagram. Similarly, chaotic firing occurs when red region in Fig. 24b and yellow region in Fig. 24c appear at the same location in the diagram.
When yellow region appears in Fig. 24b and c, the quasi-periodic firing occurs. Meanwhile, when the yellow region in Fig. 24b and blue region in Fig. 24c occur at the same location, periodic attractor will be generated.

Analysis of multiple transient transition behaviors
In Sect. 3, we have analyzed the multistability and hyperchaotic firing behaviors of Model (17)   When the initial value is selected as 1; 0; 1; À1; 1; 1 ð Þ , the System (17) shows a threestate transition, that is, the orbit of the system evolves over time with a transition of three states. Figure 27a shows the time-domain waveform of variable u when r ¼ 0:7, q ¼ 0:82, in which the two-scroll chaotic firing is located in the time interval t 2 0; 620 ð Þ, as shown in the yellow domain of Fig. 27b, c. The period-2 firing of magenta domain is located at the time interval t 2 620; 702 ð Þ . Finally, the trajectory of the System (17) tends to the blue steady-state period-2 firing. When r ¼ 0:7, q ¼ 3:1, the two-scroll chaotic  ð Þ , it evolves into another regular chaotic firing, which is shown in the green domain of Fig. 28a. Finally, the firing behavior tends to period-2 of purple domain. Corresponding phase diagrams are shown in Fig. 28b, c.

C. Four transition states
Choose the initial value is 4; 2; 1; À1; 1; 1 ð Þ , r ¼ 0:7, q ¼ 1:08; when t 2 0; 154 ð Þ, the timedomain wave presents a yellow two-scroll chaotic firing shown in Fig. 29a. The corresponding phase diagrams are shown in Fig. 29b, c. Then, the timedomain wave displays another two-scroll chaotic firing in the magenta domain when t 2 154; 185 ð Þ . When t 2 185; 251 ð Þ , the time-domain wave shows a quasi-period-2 firing in the cyan domain. Finally, the trajectory of the System (17) tends to the blue-period-4 firing behavior. When r ¼ 0:7, q ¼ 5:3, t 2 0; 122 ð Þ shown in Fig. 30a, the time-domain wave is green twoscroll chaotic firing; when t 2 122; 156 ð Þ , the time-domain wave shows a quasi-period-4 firing in the blue domain; when t 2 156; 522 ð Þ , the timedomain wave shows another two-scroll chaotic firing in the purple domain; when t 2 522; 800 ð Þ , the time-domain wave is red-period-2 firing. The corresponding phase diagrams are shown in Fig.  30b, c. To sum it up, Sect. 3.3.5 mainly introduces the transient transition behavior of the FO memristorcoupled tabu learning two-neuron model from initial state to steady state under different initial conditions when the strength of external electromagnetic stimulus changes. The evolution process of two, three, and four transition states from initial state to steady state is introduced in detail through phase diagrams and time-domain waveforms. It is worth affirming that this section points out for the first time the existence of multiple transient transition behaviors of FO memristor-coupled tabu learning two-neuron model under external electromagnetic stimulus. These studies will bring many new ideas to the theory and application of memristor in the tabu learning two-neuron model.

Application in PRNG
As we all know, chaotic systems with complex dynamics have been widely used in PRNG because of their initial sensitivity and unpredictability [50,51]. Based on the proposed FO memristor-coupled tabu learning two-neuron model, a new PRNG is designed. Other parameters remain unchanged; we generate binary sequences under the EMR strength q ¼ 5:28, initial value À0:6; À0:9; 1; À1; 1; 1 ð Þ , and the EMR strength q ¼ 1:1, initial value À1; 0; À1; 1; À1; À1 ð Þ , respectively. Binary sequences are then obtained by the following rules: (1) Given the initial conditions of the system, step size h ¼ 0:01; (2) The Adomain decomposition method (ADM) is used to iterate the FO memristor-coupled tabu learning two-neuron model at the same time, in binary form and save them as a one-dimensional array, that is, the chaotic pseudo-random sequence generated by the system.
To better evaluate the performance of the system for generating random numbers, we use the NIST 800.22 to detect and compare the randomness of twoscroll chaotic and hyperchaotic attractors under two sets of parameters and initial values. If the P-Value of NIST 800.22 is greater than 0.01, it indicates that the P-Value is evenly distributed and the sequence is random. The test report is shown in Table 1. All P-Values were greater than 0.01, so the PRNG successfully passed the NIST test and showed good randomness. Furthermore, two-scroll hyperchaotic attractor has better randomness and robustness than two-scroll chaotic attractor, that is, the high-quality  Original voice data, encrypted voice data and decrypted voice data random sequence generated by the double scroll hyperchaotic attractor is more conducive to voice encryption.

Application in voice encryption
The key length and key space of encryption method based on chaos are relatively large, especially in FO chaotic system [52]. Based on r ¼ 0:7, q ¼ 1:1 and the initial value À1; 0; À1; 1; À1; À1 ð Þ , the system generates a two-scroll hyperchaotic attractor and completes the process of voice encryption in this section. The complete block diagram of voice encryption scheme is shown in Fig. 31. The state variables S 1 ; S 2 ; S 3 ; S 4 ; S 5 ; S 6 g f of System (17) solved by ADM algorithm are normalized as follows to obtain Z 1 ; Z 2 ; Z 3 ; Z 4 ; Z 5 ; Z 6 g f : Then, the two adjacent sequences are XOR in turn to obtain five random sequences Z 1 È Z 2 ; Z 2 È Z 3 ; Z 3 È Z 4 ; Z 4 È Z 5 ; Z 5 È Z 6 g f and convert the random sequences into 32-bit binary form. Voice data are also converted to 32-bit binary form. The binary voice data and random sequence are XOR operated, and then, the XOR binary numbers are converted into floating-point numbers to generate an encrypted voice file. The decryption process is the reverse operation of the encryption process, that is, convert the random sequence and the encrypt voice sequence into binary again, decrypt it through XOR operation, and then convert the binary numbers obtained by XOR operation into floating-point numbers to obtain the decrypted original voice file.
The original, encrypted, and decrypted voice plots are shown in Fig. 32. Comparing Fig. 32a and b, it can be seen that the encrypted voice plot is very different from the original voice. Similarly, comparing Fig. 32a and c, it can be observed that the decrypted voice plot is the same as the original voice. Further analyze the error between the decrypted voice and the original voice, as shown in Fig. 33; the error is always 0, which indicates that the encryption and decryption process is successful. In addition, Fig. 34 shows the power spectrum analysis results of the encryption process. It can be seen that the power spectrum of decrypted voice is exactly the same as that of original voice. Therefore, the proposed system can be successfully used in encryption applications. In this paper, a FO memristor-coupled tabu learning two-neuron model was established by unidirectional coupling two adjacent tabu learning neurons through memristor synapse. The dynamical behaviors of the neuron model under different stimuli are further studied, that is, no external stimulus, external forced current, stimulus and EMR. In the absence of external stimuli, there are periodic attractors and transient chaos in the neuronal model. Applying external stimulus to one of the neurons can generate abundant two-scroll chaotic attractors; especially, when the neuron model is stimulated by EMR, the neuron model has multistability, hyperchaotic firing, and multitransient transition behaviors. The effects of coupling strength and external stimulation on the firing behaviors of neuronal model were studied by numerical simulation methods such as phase diagrams, Lyapunov exponent diagrams, two-parameter bifurcation diagrams, and local attraction basins. Finally, PRNG and voice encryption method based on FO memristorcoupled tabu learning two-neuron model are proposed, which have good randomness and encryption effect. In the future, we will explore more interesting memristive neural circuits and study the effects of different kinds of external stimuli with different distribution positions on the chaotic dynamical behaviors of neural networks. At the same time, according to the requirements of practical application, we will pay more attention to the physical implementation of the model.

Funding
The authors have not disclosed any funding.
Data availability Data will be made available on reasonable request.