A 6D Fractional-Order Memristive Hopfield Neural Network and its Application in Image Encryption

This paper proposes a new memristor model and uses pinched hysteresis loops (PHL) to prove the memristor characteristics of the model. Then, a new 6D fractional-order memristive Hopfield neural network (6D-FMHNN) is presented by using this memristor to simulate the induced current, and the bifurcation characteristics and coexistence attractor characteristics of fractional memristor Hopfield neural network is studied. Because this 6D-FMHNN has chaotic characteristics, we also use this 6D-FMHNN to generate a random number and apply it to the field of image encryption. We make a series of analysis on the randomness of random numbers and the security of image encryption, and prove that the encryption algorithm using this 6D-FMHNN is safe and sensitive to the key.


INTRODUCTION
The artificial neural network (ANN) system is composed of many neurons with adjustable connection weights, and has distributed information storage, good self-organization and selflearning ability, and large-scale parallel processing. ANN realizes intelligence by imitating various working and learning modes of the human brain and biological system, so it has a wide application prospect in many fields such as signal processing and pattern recognition, intelligent control, fault diagnosis, and information optimization [1][2][3][4][5].
In recent decades, more and more attention has been paid to the study of chaotic systems, and various chaotic systems with complex dynamics have been proposed [6][7][8][9][10][11]. As a new science, chaos has been widely used in secure communication [12][13][14], electronic circuit [15][16][17][18], random number generator [19][20][21][22][23], and image encryption [24][25][26][27][28][29][30][31][32][33][34]. Through the study of the brain, scientists found that there are complex chaotic phenomena in the human brain [35], which explain the irregular movement in the human brain. Inspired by this, Aihara et al. proposed a chaotic neural network model based on previous derivation and animal experiments in 1990. The chaotic neural network is a highly nonlinear dynamic system, so the neural network is closely related to chaos. Therefore, the chaotic neural network is considered to be one of the intelligent information processing systems that can realize its real-world computing. It has a good application prospect in the fields of algorithm optimization, information storage, associative memory, pattern recognition, and so on [36][37][38][39][40][41].
In 1971, Professor Leon O. Chua proposed the existence of memristor [42]. It wasn't until 2008 that HP LABS implemented the first true memristor [43]. With the physical realization of memristor, the research of memristor has attracted extensive attention again [44][45][46][47][48]. Adding memristor in the construction of nonlinear system makes the chaotic system based on memristor have richer dynamic behavior, and the generated chaotic signal has better pseudorandom characteristics, which makes it have higher research value in the fields of image encryption, spread spectrum communication, and secure communication [49][50][51][52][53]. Because memristor has memory function and its resistance value is accurate and continuously adjustable, which is similar to the function of biological neural synapse, the neural network circuit based on memristor is closer to the brain in function and can fully realize human consciousness in theory, making memristive neural network (MNN) become a new research hotspot [54][55][56][57][58][59]. In [55], the authors studied the dynamic characteristics of a small neural network with three neurons under electromagnetic radiation. The intensity of electromagnetic radiation will change the number of equilibrium points in the neural network, resulting in the diversity of attractor trajectories. In [58], a full circuit of MNN based on weighted and simultaneous disturbance training was proposed. A synaptic circuit is designed using a pair of memristors, and then a complete neural network circuit is designed.
Since the value of the next state of fractional-order chaotic system is related not only to the current state but also to all previous states, its dynamic characteristics are more complex than the integer order chaotic system. In recent years, various fractional-order chaotic systems have been widely proposed [17,[60][61][62][63]. With the deepening of research, researchers add fractional-order to the neural network model and find that the results are more similar to the activities of human brain neurons, which extends the application field of fractional-order neural network system. However, for the study of MNN, most of the results are in the integer order range [55][56][57][58][59]64]. In [65], a fractional-order memristive Hindmarsh-Rose (H-R) model was proposed. Without small parameter constraints, the fractional-order memristive H-R model had periodic and chaotic bursts, which showed that small parameters make the membrane potential activity simpler. In [66], a simplified fractional-order neural network with discontinuous conductivity function based on memristor was proposed. It was essentially a switching system with irregular switching law, which was composed of eight fractional-order neural network subsystems.
In this paper, a new memristor model is proposed and applied to the fractional-order Hopfield neural network (FHNN) to study its dynamic characteristics such as bifurcation and coexistence. In addition, the newly proposed fractional memristor Hopfield neural network (FMHNN) is also applied in the field of image encryption.

Mathematical model
In studying the relationship between four circuit variables which are voltage, current, magnetic flux, and electric charge, we know that resistance represents the relationship between voltage and current, capacitance represents the relationship between charge and voltage, and inductance represents the relationship between magnetic flux and current. According to the duality property, there should also be a circuit element that represents the relationship between charge and magnetic flux. In 1971, the concept of a memristor was proposed. This component establishes a connection between charge and magnetic flux and the relationship is described in the form of Eq. 1 [42].
Since then, many scholars have devoted themselves to the research of memristor. The concept of memristor has been extended to any device with two ends representing pinched hysteresis loops (PHL). For the first-order general memristor, its mathematical model can be described by Eq. 2 [56]: where v m and i m represent the input voltage and input current of the memristor, respectively. G is a function of x, known as memductance, which is an internal state variable. g (x, v m ) is a Lipschitz function. According to Eq. 2, we take a − b|x| − c sin(x) as memductance to construct a new memristor model Eq. 3: where a, b,and c are memristor parameters. In this paper, a, b, and c are 0.3, −0.18, and 0.2 respectively.

Pinched hysteresis loops
(2) As the frequency of the periodic driving signal increases, PHL should show shrinkage, that is, the area of PHL should constantly decrease. PHL should shrink to single-function as the frequency of the drive signal approaches infinity.
Based on Eq. 3, we draw the voltage-current relationship image of the memristor model ( Figure 1A, Figure 1B) and the time sequence diagram of voltage and current Figure 1C by MATLAB R2019 simulation software. Figure 1A shows that when the excitation voltage is v(t) sin(Ft) and the initial value is x(0) 0, the memristor exhibits PHL, and the area of PHL shrinks continuously with the increase of frequency f, which meets the two conditions of the memristor. Figure 1B shows that when the excitation voltage is v(t) A sin(t) and the initial value is x(0) 0, the PHL area of the memristor increases with the increase of the amplitude A, which is also consistent with the characteristics of the memristor.

Integer order mathematical model
Hopfield neural network (HNN) is particularly suitable for simulating various complex dynamical behaviors in the brain, especially chaotic behaviors, due to its significant nonlinearity and flexible mathematical expressions. Its mathematical model can be described by Eq. 4 [68,69]: where C i , R i , x i are membrane capacitance, membrane resistance and membrane potential of neuron i respectively. w ij is the synaptic weight between neuron i and the neuron j and tanh (x i ) is the neuron activation function. I iext is an electrical current generated by an external stimulus. According to the neural network model described by Eq. 4, HNN of three neurons can be concretized by Eq. 5 In this paper, the coefficients C 1 , C 2 , C 3 , R 1 , R 2 , R 3 are set as 1, respectively, and the weight w ij can be expressed by the weight matrix Eq. 6. It should be noted that the weight used by HNN in this paper is different from that in literature [69].
Because HNN is multiple neurons connected to each other, electromagnetic induction current will appear between HNN neurons due to the existence of membrane potential difference between neurons. In order to establish a mathematical model and study the influence of induced current between neurons on the dynamic behavior of HNN, this paper uses a magnetically controlled memristor to connect multiple neurons and uses the current flowing through the memristor to simulate the induced current between neurons [69]. Based on the above analysis, we can establish an improved HNN model. The mathematical model of Eq. 7 describes the situation in which induced current exists between three neurons in HNN model, and the bidirectional induced current is generated by simulating electromagnetic induction by interconnecting neurons with three memristors. With the addition of three memristors, a 6D memristive Hopfield neural network (6D-MHNN) model is finally established. In order to intuitively understand the 6D-MHNN model given by Eq. 7, Figure 2 presents the general schematic diagram of the extension connection between neuron and memristor in this model.

Fractional mathematical model
The concept of fractional calculus has been put forward for a long time and has become a powerful tool to study fractional differential equations and fractal functions. Since fractional order is more accurate than integer order in simulating real dynamic system, this paper will construct a fractional Hopfield neural network based on Eq. 7. At present, there are many definitions of fractional calculus, but this paper adopts Caputo's definition of fractional derivative [70,71].

Definition 1.
Caputo fractional derivative is defined as follows: where Γ(•) represents Gamma function. According to the definition of Caputo fractional derivative, the mathematical model Eq. 9 of fractional memristor Hopfield neural network (FMHNN) can be obtained by setting the order of Eq. 7 derivative as α.

System stability analysis
It is necessary to analyze the stability of equilibrium point of HNN dynamic system, so we make the left side of Eq. 9 all equal to 0, thus obtaining Eq. 10, and then solve the roots of the equations to obtain the equilibrium point.
Because it is difficult to directly solve the root of Eq. 12, the graphic method is adopted in this paper. MATLAB is used to draw the trajectory diagram of the two equations ( Figure 3A, Figure 3B, Figure 3C), and then numerical analysis is used to obtain the analytical solution of the equations.
The approximate range of x 1 , x 3 can be easily obtained from the figure, and a more accurate solution can be obtained by using Fsolve function of MATLAB. After solving the root [x 1 , x 2 , x 3 , x 4 , x 5 , x 6 ] of Eq. 10, put each solution into the Jacobi matrix of Eq. 9, and then obtain each eigenvalue from the characteristic polynomial; the results are shown in Table 1.

DYNAMICS ANALYSIS OF 6D-FMHNN
At present, related scholars usually use the phase diagram method, Lyapunov exponent method, and bifurcation diagram method to analyze the dynamics of chaotic system, etc. In addition, the attractor basin can easily find the coexistence attractor of the system. Therefore, this paper will use the FDE12 algorithm in MATLAB to solve the fractional differential equation, Eq. 9, and thus draw the bifurcating diagram, phase diagram, and attraction basin.

Dynamics of variable order
When the initial value of the system Eq. 9 is unified as [0.1, 0.1, 0.1, 0.1, 0.1, 0.1], bifurcation diagrams on the memristor coupling strength k are drawn respectively with order α = 0.7, α = 0.8, α = 0.9, α = 1 displayed on Figure 4A, Figure 4B, Figure 4C, and Figure 4D. It can be seen from the figure that when the order α is 0.7, no matter how k changes, the system will not appear as chaos.
When the order is 0.8-0.9, the value of the memristor intensity k determines whether the system behaves periodic or chaotic. When the order is 0.8, 0.07 ≤ k ≤ 0.1, the system is chaotic. When the order is 0.9, 0.007 ≤ k ≤ 0.025, or 0.07 ≤ k ≤ 0.08, the system is chaotic. When the order is one, k has to be a very small number of values for the system to be chaotic. In addition, this paper also studies the bifurcation diagram of 6D-FMHNN ( Figure 4E) on order α when the memristor coupling strength k is 0.02 and the initial value of the system is [0.1, 0.1, 0.1, 0.1, 0.1, 0.1]. As can be seen from the observation in the picture that with the increase of order, the system changes from periodic one-limit cycle to periodic two-limit cycle and then to periodic four limit cycle. Until α = 0.88, the system becomes chaotic. However, as α continues to rise, the system appears as inverse period-doubling and degenerates into a periodic four-limit cycle. After degenerating into a period, it becomes chaotic again after a = 0.93.
In conclusion, the order has a great influence on the dynamic behavior of 6D-FMHNN. When α > 0.8, the chaos phenomenon of the system is more significant.

Dynamics of fixed order
In this part, the order α of differentiation in Eq. 9 is fixed as 0.88. First, set the initial value of the system as [0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1], and explore the change of system bifurcation behavior when the memristor coupling strength k changes ( Figure 5A). Secondly, k is fixed at 0.015 and the initial value of the system was set as [x 1 (0), 0, 0, 0, 0, 0, 0], and the bifurcation diagram of the initial value x 1 of the system was drawn in Figure 5B. Observation can be found from Figure 5A that, like the general partial illustrations, with the rise of k, the system transitions from the initial period-1 limit cycle to period-2 limit cycle and then to period-4 limit cycle. When k = 0.02, the system starts to become chaotic. However, as k continues to rise, the system gradually degenerates into a periodic limit cycle of four times and bifurcates inversely when k = 0.08 and then enters a chaos state. As can be seen from the observation in the picture from Figure 5B that when the initial value of the system x 1 changes between [ − 1, 0), the system is always a periodic fourlimit cycle; when x 1 changes between (0, 1], the system presents chaos. In general, the system is sensitive to initial values and memristor coupling strength k.

A variety of coexisting attractors
Coexisting attractor means that different initial states of a system correspond to two different attractors. When the order of the system is determined to be 0.88 and the memristive coupling strength k is −0.05, the initial values of the system are [ − 0.1, 0, 0,  Figure 6D). In addition, this paper also draws the attraction basin corresponding to the initial state of [x 1 (0), x 2 (0), 0, 0, 0, 0] when the system order is determined to be 0.88 and the memristor resistance coupling strength k is 0.02 Figure 6E

Random sequence generation
Suppose the size of the image to be encrypted is M × N, then we need to generate a random sequence of length 3 × M × N, which is FIGURE 4 | The bifurcation diagram of 6D-FMHNN with respect to parameter k (A) when the differential order is 0.7; (B) when the differential order is 0.8; (C) when the differential order is 0.9; (D) when the differential order is 1 and the bifurcation diagram of 6D-FMHNN with respect to order when the memristor coupling strength k is 0.02 (E).  respectively used in the forward diffusion and reverse diffusion of the two diffusion algorithms and the scrambling algorithm.
Step 1: Let the input initial value of 6D-FMHNN system be [x 1 , x 2 , x 3 , x 4 , x 5 , x 6 ], and first iterate T wheel to eliminate the initial disturbance of the system.
Step 2: After iterating round T to eliminate the initial disturbance of the system, we conduct a slight disturbance of the initial value of every B iteration to ensure the pseudo-randomness of the sequence.
Step 3: After iterating round C, of course, we need to guarantee C × B Step 4: We reduce [X1,X2,X3,X4,X5,X6] T to a one-dimensional vector X, and then we take X(1 : 3 × M × N) so that we get the random sequence we need.

Image encryption steps
The steps of image encryption can be divided into four steps.
Step 1: By iterating the 6D-FMHNN system, the sequence with 3 × M × N length was denoted as X, and the dimension of the image to be encrypted with M × N size was reduced to a onedimensional vector.
Step 2: Make two additive mode diffusions, so that the information of each pixel point spreads to the whole image. The diffusion algorithm is shown in Eq. 13: where A is the image to be diffused, S is the random sequence with length M × N, and C is the encrypted image; so the inverse transformation of Eq. 13 is Eq. 14.
In this paper, S = X (1: M × N) is taken for the first diffusion and S = X (M × N + 1: M × N × 2) for the second diffusion.
Step 3: Shuffling pixel positions, in other words, exchanging information between two pixels via Eq. 15 mapping. Step3Perform two more multiply mode taking diffusion operations [72], and the diffusion algorithm is shown in Eq. 16: Therefore, the inverse transformation corresponding to Eq. 16 is shown below as Eq. 17: (17) where A is the image to be diffused, S is the random sequence with length M × N, and C is the encrypted image. In this paper, S = X (M × N + 1: M × N × 2) is taken for the first diffusion and S = X (M × N × 2 + 1: M × N × 3) for the second diffusion.

Image decryption step
Since the decryption operation is the reverse operation of the encryption operation, this article will not be detailed to save space.
Step 1: iterate the 6-day FMHNN system to get the sequence with length 3 × M × N, denoted as X, and reduce the dimension of the image to be decrypted with size M × N to a one-dimensional vector.
Step 2: Carry out the inverse operation Eq. 17 of multiplication for mold diffusion twice.
Step 3: Perform the same pixel scrambling Eq. 15 again.
Step 4: Perform the inverse operation Eq. 14 of two additive mode diffusions.

Randomness test of chaotic sequences
To test the randomness of the random sequence generated in part 5.1 of this article, we need to use NIST SP800-22. In this paper, the parameter k of Eq. 9 is 0.02, the order α is 0.88, and the initial value is [1][2][3][4][5][6]. NIST SP800-22 has a total of 15 test methods to test the randomness of a sequence, so the software has 15 standards for the randomness of the sequence. After the execution of each test method, a p-value between [0, 1] will be generated, and if the p-value is greater than the preset threshold, it indicates that the test has passed. In this paper, the default threshold of the software is 0.01, and the sequence length 10 6 is set to one input stream, and there are 10 such input streams. Table 2 shows the results of 15 tests.

Experimental results of image encryption
The encryption pictures in this paper come from USC-SIPI database, and the operating system of the test platform is Windows 10 64-bit, the test software is MATLAB 2019, and the processor of the test machine is quad-core Intel Core I5-7300HQ@2.50Ghz. The memory is 16 GB DDR4 RAM. The parameter k of 6D-FMHNN system is 0.02, the order is 0.88, and the initial value is [1][2][3][4][5][6]. The original, encrypted, and decrypted images of some images are shown in Figure 7. The attacker is unable to obtain any useful information from the encrypted images shown in Figure 7, indicating that the encryption algorithm proposed in this paper is effective.

Image histogram analysis
Histogram is a graph drawn after the frequency of each gray value in the image is counted. In order not to leave an attacker with any useful statistics, an ideal encryption algorithm would distribute the pixel depth of the image evenly, as shown in the histogram, so that the squares are basically flush. Figure 8 shows the histogram of pixel depth distribution of the original image and the encrypted image of some images. It can be seen from the figure that the histogram of plaintext image shows obvious statistical law, and the uniformity of pixel depth distribution of the image after encryption is significantly improved. In order to quantitatively analyze the uniformity of histogram distribution, chi-square test (Eq. 18) and variance were used to evaluate. chi-Square values of the histogram before and after encryption were calculated, and the pictures involved in the calculation include 4. . It is not difficult to see from the table that the gray value distribution of the encrypted image is quite uniform compared with that of the unencrypted image, so it is difficult for the attacker to attack by using histogram. Secondly, by changing the initial value of the encryption key, namely 6D-FMHNN system, the variance of the image histogram after encryption is calculated to evaluate whether the quality of the key has an impact on the encryption effect. It can be seen from Table 3 that the encryption  6D FMHNN and its Application algorithm in this paper has good encryption effect for different keys.
In this paper, O represents the observed count of each gray value, and E represents the expected count of each gray value.

Correlation analysis
Correlation analysis is the analysis of multiple factors with correlation, so as to measure the correlation of these factors. In this paper, it specifically refers to analyzing the correlation of the color depth of two adjacent pixels in an image. There is a strong correlation between the color depths of two adjacent pixels of an unencrypted image, and an ideal encryption algorithm should be able to hide this correlation; otherwise, it would be a breach for an attacker. In this paper, the correlation coefficient is calculated to measure the correlation of images, and the calculation formula is shown in Eq. 19 [24]: where y is the adjacent pixel of x, and N is the number of pixel pairs randomly selected from the image, that is, the correlation coefficient of the image, which is closer to 0, the better. In this paper, 5,000 pixels are randomly selected from the image and the correlation coefficients in the vertical, horizontal, and diagonal directions are calculated, respectively. The results are shown in Figure 9, Table 4, and Table 5. It can be found from Figure 9 and Table 4 that the adjacent pixels of the unencrypted picture have a strong correlation, which basically disappears after  encryption. Table 5 compares the correlation coefficients between this paper and other literatures. It can be found that our encryption algorithm has lower correlation coefficients in horizontal, vertical, and diagonal directions.

Information entropy analysis
Shannon defined information entropy as the occurrence probability of discrete random events, which is an important indicator to measure the randomness of information. It is generally believed that the greater the information entropy is, the greater the uncertainty is, so it can be used to measure the randomness of the image, and its calculation formula is given by Eq. 20 [74]: where i is the gray value, P(i) is the probability of gray value i, and L is the maximum gray value of pixel. For the image with 256 Gray level, the ideal value of H is 8. The larger the value H is, the more random the image information distribution is, and the better the encryption effect is. In this paper, the information entropy of some pictures that include boat, house, 4.  [24], [25], [73], and [74]. We can see from the comparisons that the image encryption algorithm in this paper has information entropy which is closer to the ideal value compared with other literatures.

Differential attack analysis
Differential attack means that the attacker will make subtle changes to the plaintext image, such as changing the gray value of a pixel, and then encrypt the image before and after the change, respectively, and find the relationship between the plaintext and ciphertext by comparing the two encrypted images. Therefore, an ideal encryption algorithm must be able to resist differential attacks, that is, even a change of a pixel will cause a huge difference between the encrypted images before and after. In order to quantitatively measure the ability of encryption algorithm to resist differential attack, the number of pixels change rate (NPCR) and the unified average changing intensity (UACI) are used in this paper [75]. NPCR reflects the ratio of the number of different pixel values in the same position of two images to all pixels, and UACI reflects the average change intensity of the difference between the data of two encrypted images. The calculation formula of the two indicators is shown in Eq. 21 [25]: where M and N are the size of the image, T is the maximum gray value of the image, and C 1 C 2 is the original encrypted image and slightly modified encrypted image, respectively. The ideal values for NPCR and UACI are 99.609% and 33.463 5%, respectively. Using Eq. 21, UACI and NPCR of a series of images were calculated in this paper and compared with other literatures. The results are 99:6075% and 33:4667%, respectively. In [24], the results are 99:6352% and 33:5,614%. In [25], the results are 99: 609% and 33:4907%. In [73], the results are 99:6114% and 33: 4523%. In [74], the results are 99:5,956% and 33:4535%. It is obvious that the image encryption algorithm in this paper has UACI/NPCR which is closer to the ideal value compared with other literatures.

Key sensitivity analysis
The general method for attackers to crack keys is brute force cracking, that is, to try all possible keys one by one. Because existing encryption algorithms generally have sufficient key space, which means that brute force cracking takes a long time; even for computers with powerful computing power, brute force cracking is not practical. However, if the encryption algorithm does not have key sensitivity, in other words, the ciphertext obtained by encrypting the same image with two similar keys is also similar, the attacker will take advantage of this vulnerability, optimize the brute force cracking algorithm, and quickly crack the key. In this paper, key1 = [1-6] and key2 = [1 + 10 -10 , 2, 3, 4, 5, 6] are used to encrypt the same image, respectively, to obtain ciphertext 1 and ciphertext 2 (Fig. 12), and the NPCR and UACI values of ciphertext 1 and ciphertext 2 are calculated as NPCR = 99.612 4% and UACI = 33.441 0%, respectively, which means that the two ciphertexts have great differences. In addition, key2 is used in this paper to decrypt the ciphertext of key1, and the results obtained are shown in Figure 10. Therefore, this algorithm has good key sensitivity.

CONCLUSION
This paper proposes a new memristor model and uses PHL to prove the memristor characteristics of the model. In addition, we propose a new 6D-FMHNN by using this memristor to simulate the induced current caused by potential difference between neurons and studies the dynamic behavior of 6D-FMHNN, such as bifurcation characteristics and coexistence attractor characteristics. 6D-FMHNN is chaotic and periodic due to different order, parameter, and initial value. Because this 6D-FMHNN has chaotic characteristics, it can be used to generate random sequences, so we also use this 6D-FMHNN to generate random numbers and apply them to the field of image encryption. This paper makes a series of analysis on the randomness of random numbers and the security of image encryption, and proves that the encryption algorithm using this 6D-FMHNN is safe and sensitive to the key. In the future, we will look for FMHNN that can generate multiple scrolls or multiple attractors.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
FY, XK, and SC contributed to conception and design of the study. HC organized the database. QY performed the statistical analysis. XK wrote the first draft of the manuscript. FY, XK, SC, and SD wrote sections of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.