Generation of Multi-Lobe Chua Corsage Memristor and Its Neural Oscillation

The Chua corsage memristor (CCM) is considered as one of the candidates for the realization of biological neuron models due to its rich neuromorphic behaviors. In this paper, a universal model for m-lobe CCM memristor is proposed. Moreover, a novel small-signal equivalent circuit with one capacitor is derived based on the proposed model to determine the edge of chaos and obtain the zero-pole diagrams and analyze the frequency response and oscillation mechanism of the m-lobe CCM system, which are discussed in detail. In view of existence of the edge of chaos, the frequency response and the oscillation mechanism of the simplest oscillator is analysed using the proposed model. Finally, the proposed model has exhibited some essential neural oscillation, including the stable limit cycle, supercritical Hopf bifurcation, spiking and bursting oscillation. This study also reveals a previously undiscovered behavior of bursting oscillation in a CCM system.


Introduction
The development and complexity analysis of biological emotion-/memory-like models, especially the neural oscillation and design of the related mimicking electronic units, have become the hot topics in the fields of neuromorphic engineering and neuroergonomics. Due to their nature of unique memory, the memristors and memristive devices have become one of the strong candidates to achieve, simultaneously, scalability and biological flexibility, which push forward the next generation of neuromorphic computing. Therefore, as one of the general memristors, Chua corsage memristor (CCM) has been considered the key element to imitate the nonlinear biological behaviors and construct the neural network [1]. In 2016, the small-signal equivalent circuit and the simplest electronic oscillator consisting of only one CCM in parallel with a battery was presented [2]. Subsequently, a series of studies was conducted to investigate the complex characteristics of CCMs, such as Hopf bifurcation [2], the locally active [3,4], the edge of chaos [5], pinched hysteresis loops [6,7], basin of attraction [8,9], phase portraits basin of attraction, phase portraits [10], etc.
Hereby, some conclusions in the existing research literature could be reviewed as follows: (i) the theoretical studies and dynamics analysis on 2-/4-/6-lobe CCM systems [7,[11][12][13][14] have been implemented by H. Kim and his team. Subsequently, the global dynamics, locally active, phase portraits, basin of attraction, pinched hysteresis loops, switching kinetics, physical realization, oscillation, and their emulator circuits [15,16] were analyzed. However, most of these research focused on several specific CCM systems. When the more lobes are needed, some important questions might be asked. For example, whether these typical nonlinear behaviors still occur in the other family of CCM, such as 12-or 14-lobe, even 20-lobe ones. What are the differences between commonality and individuality for different m-lobe CCMs? Since the small-signal equivalent circuit and the admittance function with one inductor has been designed in 2015, has other types of equivalent circuits with one capacitor and their impedance functions been developed and considered? Whether more than one similar oscillator could be built; (ii) the neural oscillation phenomenon have been found in some memristors and memristive systems, such as memristive synaptics, memristor-based neurons [17][18][19][20] and neural networks [21][22][23][24].
Then, the following studies should be worth pondering: whether most of existing models could be constructed based on the m-lobe CCM? When the biological neurons and neural networks needed to be simulated, why the m-lobe can be thought one of the powerful candidates? (iii) the applications, design, and physical implement of the CCM have also been carried out, such as the memristive self-learning logic circuit [25], in-memory computing [26], bistable nonvolatile [27], tri-stable locally-active, neuromorphic dynamics [28], neural oscillation dynamics [29,30], etc., but still in their infancy.
Based on the above literature review, the following topics with the aims of furthering research on the CCM are presented in this paper: (i) the universal model of the m-lobe CCM is introduced, which is helpful to answer the differences between commonality and individuality for all the families of CCMs; (ii) in order to investigate the versatility and flexibility for the proposed model, universal impedance function, and a novel small-signal equivalent circuit with one capacitor is designed, which is helpful for understanding the electronic circuit theory and analyzing dynamics behaviors; (iii) based on the existing theory, the nonlinear dynamics and neural oscillation on the proposed model are illustrated, such as the edge chaos, neural oscillation (i.e., supercritical Hopf bifurcation, spiking and bursting oscillation.), nonvolatile, coexisting pinched hysteresis loops, etc.
The remainder of this paper is organized as follows: in Section 2, the universal model of m-lobe Chua corsage memristor is introduced. We also summarize several natural characteristics of the proposed model, such as pinched hysteresis loops, DC V − I curves, and multiple locally-active domains. In Section 3, the universal impedance function is derived. Then, the novel small-signal equivalent circuit with a single capacitor is designed. The existence of the edge of chaos, zero-pole diagrams, and frequency response are discussed and summarized. Moreover, in Section 4, the distributed rules on the external positive inductance are desmonstrated. The oscillation mechanism, supercritical Hopf bifurcation, limit cycle, spiking and bursting oscillation of the simplest oscillator circuit is demonstrated based on the relationship between the model parameters and the excitation voltage. Finally, the paper is summarized in Section 5.

The Universal Model of m-Lobe Chua Corsage Memristor
The universal model of m-lobe CCM can be introduced as follows: G(x) is the memductance function and G 0 (G 0 = 0.1) is a real constant; i, v and x are denoted as the current, voltage and state variable of the CCM, respectively. A and a = (N − 1) 2 + 1, (N is an integer, N = 2, 3, 4, . . . ) are the real parameters; A M is called the midpoint briefly; both n j and h j are the appropriate integers (h j > n j ); m = 2, 4, 6, . . . represents the number of the lobe. Then, the coexisting pinched hysteresis loops are shown in Figure 1.

Proposition 1.
For the proposed universal model, the parameter A exhibits three cases, which imply the quite different conditions and negative slope domains (namely, locally active regions), the related statement as follows: Case 1: When the parameter A = A M is fixed, only one negative slope domain exists in the first stable branch, which is located over the negative range of the voltage; Case 2: When the parameter A = A FS < A M is chosen, only one negative slope domain could be observed in the first branch, which is located over the positive range of the voltage; Case 3: When the parameter A = A FL > A M is obtained, two the negative slope domains could be obtained in the first and second stable branch.

Proposition 2.
Several natural characteristics of the proposed models in both v vs. i and dx/dt vs.
x are shown as follows: (1) m + 1 equilibrium (i.e., Q 0 , Q 1 , Q 2 , . . . , Q m ) and m turning points (namely., T 1 , T 2 , T 3 , . . . , T m ) could be observed in their dynamic route maps (DRMs); (2) Substitute the x-axis coordinate of the first DC operating point (T 1 ) into f (x), then f (T 1 ) < 0; (3) Substitute the x-axis coordinate of the last DC operating point (T m ) into f (x), then f (T m ) > 0; (4) The properties and nonlinear phenomena (i.e., the location of DC operating points, local stability, active and locally active regions, oscillation, etc.) would not disappear as the number of lobes increase; (5) Edge of chaos can occur in the first lobe, but do not change and disappear as the number of lobes increase; (6) Similar neural oscillation dynamical behaviors might emerge. However, they are not identical phenomena in all the m-lobe CCMs.
Notably, before analyzing the proposed universal CCM model, some statements need to be elaborated: (i) the above features (4)-(6) would be verified via Figures 2-7. Observed from these Figures 3-9, the negative slopes are mainly concentrated in the first lobe; (ii) there are plenty of ways to set the values of the parameters n j and h j . In this paper, one method is chosen to obtain the m-lobe CCM, i.e., n j = (j + 1)(j + 2), h j = n j+1 , (j = 1, 2, . . .); (iii) since the complexity and dynamics oscillation are mainly concentrated in the first lobe, a 2-lobe CCM as an example can be considered an example, the verification and discussion on the above characteristics are given as follows.  When driven by one periodic input voltage/current excitation with a zero DC component, the frequency-dependent pinched hysteresis loops become a signature. Next, the Lissajous figures of the proposed model in the i − v plane with a sinusoidal input signal are exhibited in Figure 1.
As expected, from Figure 1, the pinched hysteresis loop for amplitude, A m = 1, exists for all frequencies, ω = 0.1 rad/s, 5 rad/s, 10 rad/s, 90 rad/s, and 100 rad/s. All pinched hysteresis loops pass through the origin. Besides, the lobe areas of the pinched hysteresis loops shrink as the frequency increases [31], thereby their fingerprints are confirmed.

Parametric Representation and DC V − I Curve
According to the description on the obove natural characteristics, the number and location of locally active domains are determined by the parameter A. Herein, the there cases (i.e., A = A M , A = A FS < A M , and A = A FL > A M ) are discussed and their DC V − I curves be drawn as follows: , the simplest 2−lobe CCM can be obtained and rewritten as: Then, the dynamic route maps (DRM) of the model (2) are shown in Figure 2.  From Figure 2a, when 2-lobe CCM is short-circuited (v = 0), the direction of motion of the state variable x from any initial state x(0) is indicated by the arrowheads on the poweroff-plot (POP). The features of two turning points (namely, T 1 , T 2 ) and three equilibrium points (i.e., Q 0 (x = 3), Q 1 (x = 9), Q 2 (x = 15)) are identical with the description in Proposition 2 (1)∼(3). Both Q 0 and Q 2 are stable equilibrium points, whereas Q 1 is unstable due to the state variable x(t) diverges away from Q 1 . The dynamic routes for v = 0, ±5, ±10 are depicted in Figure 2b. Each DRM route is divided into 3 segments by breakpoints at x = 6 and x = 12, which has two asymptotically stable and one unstable equilibrium points. Moreover, two enclosed areas can be observed from Figure 2a, i.e., areas ∆ 011 and ∆ 122 formed by points Q 0 , T 1 , and Q 1 and Q 1 , T 2 , Q 2 respectively.
Generally, the DC V − I curve is used to make sure the basic parameters for the proposed CCM. Therefore, the explicit formula is obtained from Equation (2), which equals zero and computes G 0 = 0.1, x = X as the functions of v = V and i = I, The corresponding DC V − I curves are drawn in Figure 3 over the input variable range x ∈ [−10, 20] and voltage range V ∈ [−10, 5].
In Figure 3, abserved from between Figure 3a,b, the DC V − I curve emerges only one negative slope domain in the first stable branch over the positive variable ranges 0 < x < 2 and the negative voltage range −3 V < v < −1.0 V, which lies in the third quadrant of the V m − I m coordinates and implies the existence of the locally active region. It satisfies the description in case (1). Besides, in Figure 3b, the slope at in blue; the slope V r = 2.9 V is 10 (got from the coordinates at V r = 2.91 V and V r = 2.89 V) in red.
Observed from Equation (1) and Figure 3, when the parameters A = 30, n j = 20, h j = 40 are chosen, the proposed universal model is equivalent to the classical CCM model, which satisfies the rule A = A M and conforms to the description on its natural features. That means their DC V − I curves have two stable branches and one unstable branch, and only one negative slope domain in the first branch over the negative voltage range; both upper left lobe and lower right lobe have similar shapes and areas (called a symmetrical "Corsage ribbon"). Hence, they could be considered as belonging to the same family of universal 2-lobe CCM.
Case 2: A = A FS < A M In this case, we choose a new group of parameters (n j = 6, h j = 22, m = 2, A = A FS < A M = 9), the proposed model is restated as follows: The DRM for the POP and v = 0, 5, 10 are plotted in Figure 4. Comparing Figure 4a with Figure 2a, the area ∆ 011 enclosed by three points (i.e., Q 0 , T 1 , and Q 1 ) is much larger than ∆ 122 by three points (i.e., Q 1 , T 2 , and Q 2 ), which implies the different shapes and areas of the corsage ribbon.
With the new parameters, the graphs of Ivs.X and DC V − I are redrawn and shown in Figure 5. There is only one negative slope domain in the first branch over the positive voltage ranges 2.3 V < v < 7.0 V in Figure 5b, which lies in the first quadrant of the V m − I m coordinates and implies the existence of locally active regions. It satisfies the description in case (2)  Observed from Figure 5, the upper left lobe becomes much bigger than the lower right one, which is called the asymmetrical "Corsage ribbon".
The third group of parameters (n j = 6, h j = 22, m = 2, and A = A FL = 20 > A M ) are given for showing the more active and locally active domains clearly, the proposed model is recast as follows: The DRM for the POP and v = 0, 5, 10 are plotted in Figure 6. Comparing Figure 6a with Figure 2a, the area ∆ 011 enclosed is much smaller than the area ∆ 122 , which is also different from the previous cases. It is not only multi-valued, very unusually, but resembles an asymmetrical phenomenon. Moreover, the graphs of the Ivs.X and DC V − I curves are depicted in Figure 7.
It can be seen from Figure 7, there are two negative slope domains in the stable branches over the ranges of −4.0 V < v < −1.33 V and −14.0 < v < −12.0 V via the red and blue solid lines, which indicate the existence of two locally active regions and both lie in the third quadrant of the V m − I m coordinates. This situation satisfies the description of the case 3 in the above section. Additionally, several important slopes can be calculated 12 Siemens in the red point. Therefore, the asymmetrical "Corsage ribbon" can emerge.
To summarize, the curves of DRMs and phase portraits with different cases (i.e., A FS , A M , and A FL ) are illustrated in Figure 8.
Based on the previous description, the following conclusion can be drawn. When the other parameters (such as m, n j , and h j ) are fixed: (1) the shape and symmetry of the lobes, stability, locally active domains depend on the value of parameter A; (2) in order to obtain the negative slope domains, the ranges of the variables x and v at the equilibrium point (Q) obey the condition: xv < 0.

The Generation of Multiply Lobes for the Universal CCMs
The proposed m-lobe CCM can be implemented by model (1). The number of lobes can be captured and their parameters are listed in Table 1.
Additionally, the generation of m-lobe CCMs are shown in Figure 9. According to the above description, it can be seen that the local active domain (LAD) in the first stable branch lobes might be impacted but will not disappear with the increase of the number of lobes. Moreover, for any m-lobe CCM, the feature of parameter A will not change as the number of lobes increases, which is in line with statements on the natural characteristics of the CCM.  = 1, 2, 3, . . . ) f (x, v) Figure   2-lobe 1 2 0 9 n 1 = 6, h 1 = 12 9 − x + |x − 6| − |x − 12| Figure 3b 4-lobe 2 4 2 17 Figure 9a 6-lobe 3 6 5 28 10-lobe 5 10 17 57 12-lobe 6 12 26 75 14-lobe 7 14 37 97

Small-Signal Equivalent Model with One Capacitor
In order to predict the response of a memristor to a small-signal input applied at an equilibrium point, the small-signal equivalent of the CCM has been designed by L.O Chua and his coworkers in 2015. They have also pointed out that adding at least one energy storage element across the Chua Corsage Memristor could make sure the circuit oscillates. In these papers, the analysis on an inductance and two resistances were elaborated. Moreover, they have verified that a positive inductance is needed to compensate for the imaginary part of admittance Y(iω, V) as well as to make the total impedance to zero. As standard electronic circuit theory, some conclusion can be directly obtained via impedance function instead of admittance, which is our focus in this subsection.
The small-signal impedance (Z(s, V)) at the equilibrium points can be summarized and uncovered as follows: where a 11 , a 12 , b 11 , b 12 express the parameters of the Laplace transformation for δi m (t) and δv m (t), which neglects the higher-order terms by |δi| 1 and |δv| 1. Then, The following equivalent circuit form is rewritten from Equation (6): Then, the small-signal equivalent circuit with one capacitor and two equivalent resistors can be designed, Besides, the schematic diagram of the small-signal equivalent circuit of the impedance function with a capacitor (C m ) and two resistors (R m and R n ) is drawn in Figure 10. Then, the curves defining the parameters of this novel small-signal circuit model are shown in Figure 11. Both C m and R m are positive and R n is negative. Since the small-signal impedance Z(s, V) at v = V has only a real zero s = b 11 (V), which is necessary to add at least one capacitor (or inductor) across the Chua Corsage Memristor in order to make the circuit oscillate at some frequency ω 0 > 0. Its universal impedance function can be expressed as follows where, k = R m R n /(R m + R n ) is represented as the gain of impedance, p = −1/[(R m + R n )C m ] is denoted the only real pole, and z = −1 is only real zero for the proposed universal m-lobe CCM.
To determine the type of energy-storage element (inductor or capacitor), the following frequency response Z(iω, V) should be derived and plotted by substituting s = iω. The frequency response of impedance function, Re[Z(iω, V)] and Im[Z(iω, V)] are given in Equation (11): Observed from Figure 11, we can clearly see that the configured small-signal equivalent circuit and the analysis on its impedance function verify the oscillation theory and the existence of oscillators from another perspective.

Edge of Chaos
Based on the representation of zeros and poles in the impedance function, the proposed model clearly exhibit the edge of chaos as shown on the zero-pole diagrams. Observed from Figures 3, 5 and 7, different negative slopes regions can be captured by adjusting parameter A. The relationship of parameters x and v conforms to the rule xv < 0.
Besides, the feature of parameter A can be summarized as: when A ≤ A M , there exists only one negative slope region; when A > A M , there are two.
Next, the gain, zero and pole of the Equation (9) are calculated and the diagrams are plotted in Figure 12 to exhibit the edge of chaos. From Figure 12, the distribution of LADs are identical with the DC V − I curves in Figures 3, 5 and 7. There is also one energy storage element that is indeed required for generating an oscillator, which could be an inductor or capacitor in serial/parallel with the resistor.

The Simplest Oscillator and Its Neural Oscillation
In light of the concepts and techniques of both electric circuit and nonlinear dynamics theory, there is a pair of complex-conjugate poles on the imaginary axis at some frequency ω > 0 in one oscillator circuit. According to [1][2][3]11,14], one external inductance (L * ) can be chosen to present plenty of nerual oscillation. In this section, the discussion on the properties and rules for this positive inductance are given as well as its frequency response and neural dynamics behaviors.
In Figure 13, L * is an external positive inductance, v M is considered the voltage of the proposed universal m-lobe CCM; V represents a DC voltage source.

The External Positive Inductance
Based on the frequency response of the derived impedance function (Z(iω, V)) and the relationship between admittance and impedance (Y(iω, V) = 1/Z(iω, V)), the external inductances (L * ) can be calculated in three cases.
Case 1: A = A M The frequency response for the model (2) are illustrated in Figure 14. The admittance Y are shown in Figure 14a over the frequency range −50 rad/s ≤ ω * ≤ 50 rad/s.

Case 2: A = A FS < A M
In this case, the frequency response of the model is illustrated in Figure 15. The complex admittance Y(iω * , V) are plotted in Figure 15a over the frequency range of −50 rad/s ≤ ω * ≤ 50 rad/s. The Re[Y(iω * )] = 0 at ω * = ±9.51 rad/s whereas Im[Y(iω * )] = ±0.1046 with the external L * are calculated as follows: Observed from Figure 15b, when L * = 1.005H is chosen and the voltage lies in the LAD, the imaginary part can be compensated. Some parts of the loci shown in Figure 15b are located in the open RHP, which implies that the circuit can be destabilized by varying the inductance L * to generate the desired nonlinear behaviors.

Case 3: A = A FL > A M
In this case, the frequency response curves obtained from both LAD1 and LAD2 are shown in Figure 16. For v in LAD1, the oscillation occurs when ω * = 7.27 rad/s with Y 1 = 0 ± 0.1356i and ω * = 1.83 rad/s with Y 2 = 0 ± 0.5463i, respectively. The corresponding inductance L * for both LADs can be calculated by the following formulas: According to Equations (13)-(15), the following conclusion can be reached: in the simplest oscillator circuit, in order to compensate for the imaginary part of admittance and capture the nonlinear behaviors, the external positive inductance (L * ) always need to be connected in series with all types of the small-signal equivalent circuits. The number and distribution of L * is determined by the parameter A and the LAD(s).
Subsequently, it is summed up as follows: when the parameter A ≤ A M , only one external inductor (L * ) is needed to compensate the imaginary part of admittance and generate the neural oscillation. Whereas, when the parameter A > A M , more than one inductance value can be used to design the oscillator.

Oscillation Mechanism
The state equations of the simplest oscillator circuit in Figure 13 can be written as follows: where A, n j , and h j denote the parameters; x and i are the state variable and the current of the inductor (L * ); L * is a positive inductance, and G 0 = 0.1. Then, the output voltage (v out ) is set as the voltage of the universal m-lobe CCM, that is, v out = −v M .
To ease the demonstration of local activity and edge of chaos, from which we can identify the mechanism of the action potential in this circuit, the following universal equation of Y 2 (s, V) at the equilibrium point Q for the oscillator circuit is formulated: where the impedance of an external positive inductance (L * ) is represented as Z L * = L * s. The zeros and pole of Y 2 (s, V) are solved as follows: From Equation (18), both poles p 21 and p 22 are the functions of the external positive inductance (L * ). When the condition (R n R m C m + L * ) 2 + 4R n L * /z 2 < 0 is satisfied, an oscillator emerges. Correspondingly, the frequency response can be moved from the open LHP to the open RHP by increasing the inductance L * ∈ [0, +∞).
According to the admittance function and frequency response of Figure 13, Hopf bifurcation, stable, and destabilized phenomena might exist, as well as the spiking and neural dynamics, which are the main focus in the next subsection.

Oscillation and Neural Dynamics
(1) Limit cycle and supercritical Hopf bifurcation.
Hopf bifurcation gives birth to a limit cycle to change the nonlinear system stability [2,13]. A stable limit cycle could lead to the supercritical Hopf bifurcation. Therefore, for the proposed universal CCM, its oscillators can exhibit the supercritical Hopf bifurcation as shown in Figure 17 for A = A M , Figure 18 for A = A FS , Figures 19 and 20  Observed from Figures 17-20, the limit cycle lies in the open RHP of its pole plot between two Hopf bifurcation points. Therefore, it can be confirmed that the stable limit cycle does appear.
(2) Neural oscillation When the memristive emulator is utilized to mimic one biological neuron, spiking oscillation can be considered one of the most prominent phenomena. During the supercritical Hopf bifurcation intervals, the large-amplitude non-sinusoidal periodic waveform can be discovered in three cases, which is shown in Figure 21. It is evident that the domain of spiking periodic oscillation could occur depending on the parameters of the proposed universal CCM model, which is not mentioned in the supercritical Hopf bifurcation theory. Finally, the non-periodic oscillation with a small amplitude could exist on the open LHP that that is far away from the Hopf bifurcation points in Figure 21d. In other words, the non-sinusoidal "spiking" periodic oscillation, periodic and non-periodic oscillation can coexist in the proposed universal model (1); (3) Bursting oscillation From above figures, the distributions of the spiking oscillation (i.e., action potential) has been analyzed, which resemblance to biologically-generated action potentials at the same time, it verifies the Chua's theory on "local activity" is the origin of action potential (spikes). In this subsection, in order to further observe the neural dynamics for the proposed CCM-based second-order circuit, the other significant nonlinear behavior is presented, such as bursting oscillation and their distributions, which are illustrated in Figure 22. Observing all the bursting regions. They lie on the open LHP with a non-periodic small-amplitude oscillation, but unstable, and rapidly transition into a stable (oscillating) motion as V increases in LAD(s). Among them, the widest region is the case of A = A FL , and the narrowest one is A = A M . There are two regions that have been found in the case of A = A FL . In other words, more complex nonlinearity and neural dynamics exist in the case of A = A FL , such as edge of chaos, local active domain, the compensated positive inductance (L * ), limit cycle, supercritical Hopf bifurcation, spiking and bursting oscillation, and so on. As the essential of the bursting which implies a series of neural behaviors, which is one of the most important topics in the research of neuromorphic, such as biological adaptive behavior, long-short-term memory, biological emotion-/memory-like models, etc. Although all the bursting oscillations are distributed in a tiny-region(s), they still could be used to construct the emulated electronic neuron or neural network circuits. Clearly, the bursting phenomenon agrees with the engineering requirement in the field of biomimetic neurology, and is our primary contribution for the proposed CCM.

Conclusions
In order to gain in-depth understanding of the nonlinear characteristics of the Chua corsage memristor and explore the applications in the field of biomimetic neurology, the universal model and generation rules of the m-lobe Chua corsage memristor are introduced, as well as their natural characteristics. The novel small-signal equivalent circuit with one positive capacitance and two resistors is proposed and its impedance function is presented and analyzed to verify the related theory. Moreover, the zero-pole diagrams and frequency response of the admittance functions for the simplest oscillator are discussed in detail, according to the existence of the edge of chaos. Furthermore, the features and distribution of the compensated positive inductance (L * ) are discussed and summarized. In addition, the oscillation mechanism of the proposed CCM-based simplest oscillator is analyzed. Finally, the neural dynamics are demonstrated, such as the limit cycle, supercritical Hopf bifurcation, spiking and bursting oscillation, respectively. This study provided a theoretical foundation for application in the field of biomimetic neurology.