Turing instability and pattern formation of a fractional Hopﬁeld reaction–diffusion neural network with transmission delay

. It is well known that integer-order neural networks with diffusion have rich spatial and temporal dynamical behaviors, including Turing pattern and Hopf bifurcation. Recently, some studies indicate that fractional calculus can depict the memory and hereditary attributes of neural networks more accurately. In this paper, we mainly investigate the Turing pattern in a delayed reaction–diffusion neural network with Caputo-type fractional derivative. In particular, we ﬁnd that this fractional neural network can form steadily spatial patterns even if its ﬁrst-derivative counterpart cannot develop any steady pattern, which implies that temporal fractional derivative contributes to pattern formation. Numerical simulations show that both fractional derivative and time delay have inﬂuence on the shape of Turing patterns.


Introduction
In the past few years, neural networks have attracted much attention in different fields of science and engineering, owing to their valuable applications in associative memory, combinatorial optimization, image processing and so on [5,13]. Some of the above applications are based on the dynamical behavior analysis of neural networks [2]. Actually, neural networks are realized by large-scale integrated circuits, and the density of electromagnetic fields is generally not uniform. Therefore, in factual modeling, only considering the change of time seems to be not comprehensive when electrons are moving in asymmetric and nonuniform electromagnetic fields [21,22]. Due to the reaction-diffusion effect, neural networks have rich spatial and temporal dynamical behaviors, like various Turing patterns.
Spatial dynamics in reaction-diffusion systems was originally proposed by Turing in the work of [23], which mainly focuses on chemical reaction systems. This pioneer work not only established the theoretical foundation for understanding diverse patterns occurring in the natural world, but also opened a new research field, namely, pattern dynamics, which has received extensive attention and is still a hot topic in many fields such as species dynamics [7,27,28], medicine [26,31], networks [4,14,29]. Based on the reaction-diffusion theory of Turing, Chua and Goraş [4] investigated the phenomenon of pattern formation in cellular neural networks. Recently, Zhao et al. [29] studied the conditions of Hopf bifurcation and Turing instability in a reaction-diffusion neural network, where the corresponding model is described as follows: ∂u(x, y, t) ∂t = d 1 ∆u(x, y, t) − c 1 u(x, y, t) + a 1 tanh v(x, y, t) + b 1 tanh u(x, y, t) , ∂v(x, y, t) ∂t = d 2 ∆v(x, y, t) − c 2 v(x, y, t) + a 2 tanh u(x, y, t) + b 2 tanh v(x, y, t) (1) in which u(x, y, t), v(x, y, t) stand for state variables of neurons at time t and spatial position (x, y); c i > 0 (i = 1, 2) denote the rates of resetting neuronic potential to the resting state in isolation; a i , b i (i = 1, 2) represent connection weights; ∆ = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 is the Laplacian operator in a two-dimensional space; d i > 0 (i = 1, 2) are diffusion coefficients of electrons among neurons. Traditional neural networks, such as system (1), are mainly established from the view of integer-order derivatives, which can be described by classical ordinary differential equations. From the perspective of the application the integer-order derivative is used to describe some properties at a certain time in physical processes or some local properties of a certain position. In recent years, experimental research indicates that fraction-order derivatives provide an excellent tool for the description of memory and hereditary properties of various materials and processes [1,15]. Generally speaking, plenty of practical objects can be described clearly by the fractional differential equations due to their more degrees of freedom and infinite memory. Thus, the research of fractional neural networks has gained much attention, and some valuable results have been referred to [9][10][11][12]20].
Besides, due to the finite speed of signal transmission and amplifiers switching, time delay exists unavoidably in neural networks, which is known as transmission delay. Usually, the time delay is harmful to the stability performance of neural networks, causing oscillation, even chaos. Zhao et al. [30] introduced the transmission delay into a discretetime Hopfield neural network model and showed that Hopf bifurcation occurs when the delay exceeds a critical value.
In [8], fractional calculus is used in electrical circuits that possess memory and hereditary properties, which is similar to the integrated circuits of Hopfield neural network.
Moreover, since fractional calculus is a generalization of conventional calculus, it is expected that fractional models will generally provide a more accurate description of the system dynamics than those based on classical differential equations [17]. In order to describe the memory and hereditary properties of Hopfield neural network, fractional order q is introduced into our model. Besides, time delay τ 1 and τ 2 are considered, which accounts for the finite speed of signal transmission and amplifiers switching in integrated circuits. Furthermore, it is generally assumed that the reaction-diffusion terms in system (1) play a conclusive role in pattern formation. A question arises naturally whether temporal fractional derivative has some relationship with pattern formation. To be specific, can the fractional reaction-diffusion neural network produce steadily spatial patterns even if its first-derivative counterpart cannot form any steady pattern? To this end, we consider the following delayed fractional reaction-diffusion neural network: with initial condition and Neumann boundary condition where q ∈ (0, 1] is the fractional order; τ 1 and τ 2 are transmission delay that reflects the time-lag effect; f (·) denotes the activation function, which maps the input to the output of the neuron; a square domain Ω = (0, L) × (0, L) in which L is a positive bounded constant; n is the outward unit normal vector of the boundary ∂Ω that is assumed to be smooth; other variables and parameters are similar to those in system (1). For simplicity, let τ = τ 1 = τ 2 . The Neumann boundary condition implies that nothing enters this system and nothing exits from this system. The main contributions of this paper can be summarized as follows: (i) By theoretical analysis and numerical simulations, the paper illustrates that temporal fractional derivative contributes to pattern formation in system (2); (ii) Fractional derivative can soften the stability conditions of the homogeneous steady state, which implies that when fractional-order q = 1, Hopf bifurcation occurs in system (2), but once q < 1, the homogeneous steady state may become stable; (iii) To avoid complicated and tedious MATLAB programming, numerical simulations are mainly completed by Simulink, which is more efficient and can visually represent the interrelationships of neuron states; (IV) Numerical simulations indicate that the fractional derivative is related to the shape of the Turing pattern. Besides, as the fractional derivative decreases, it softens the Turing instability conditions of system (2). This paper is organized as follows. In Section 2, we first analyze the existence of Hopf bifurcation and obtain the parameter sets, where the homogeneous steady state of system (2) is stable. Furthermore, we present and discuss the results of Turing instability and pattern formation, which will be illustrated by numerical simulations in Section 3. The paper ends with a conclusion in Section 4.

Hopf bifurcation and Turing instability
Throughout this paper, we address the following assumption holds.

Hopf bifurcation
The homogeneous steady state of system (2) is E 0 = (0, 0). In this subsection, we focus on the local asymptotic stability of E 0 , which is useful for Turing instability analysis. The linear perturbation equations with respect to E 0 are where λ is the growth rate of perturbations in time t, i is the imaginary unit and i 2 = −1.
In order to facilitate the analysis, let k = k 1 = k 2 . Substituting (4) into (3), we obtain the characteristic equation as follows: where To find possible periodic solutions, which may bifurcate from a Hopf bifurcation point, let λ = iω (ω > 0 is a real number) be a root of (5). Separating the real and imaginary parts yields It follows that From the first equation of (6) we have for n = 0, 1, . . . , Define the bifurcation point of system (2) as τ 0 = min{τ (n) }. To establish the main results of this section, we list the following assumptions.
Define a parameter set where E 0 is locally asymptotically stable when τ = 0 and q = 1. For τ = 0 and q ∈ (0, 1), we first introduce the following lemma.
Lemma 1. (See [16].) The following autonomous system In this case, each component of the states decays towards 0 like t −q . Also, this system is asymptotically stable if and only if | arg(λ i )| qπ/2 and those critical eigenvalues that satisfy arg(λ i )| = qπ/2 have geometric multiplicity one.
From Lemma 1 we obtain that the homogeneous steady state E 0 is locally asymptotically stable if the following inequalities hold: |arg(λ i )| > qπ/2, i = 1, 2. Similarly, we define the following parameter set: where E 0 is locally asymptotically stable for τ = 0 and q ∈ (0, 1).
In conclusion, when τ = 0, the homogeneous steady state E 0 is stable if the parameters ϕ 1 and ϕ 2 are in the set Π 1 ∩ Π 2 .
The homogeneous steady state is locally asymptotically stable for τ ∈ [0, τ 0 ). (b) System (2) undergoes a Hopf bifurcation at the homogeneous steady state when τ = τ 0 , i.e., it has a branch of periodic solutions bifurcating from the homogeneous steady state near τ = τ 0 .
Choose (ϕ 1 , ϕ 2 ) as free parameters and project the parameter sets Π 1 and Π 2 into a two-dimensional space with respect to (ϕ 1 , ϕ 2 ). From the conditions of parameter sets Π 1 and Π 2 we draw the stability regions of E 0 . In Fig. 1, the parameter set Π 1 corresponds to region A (green), while Π 2 corresponds to the region B (yellow). Region B implies that temporal fractional derivative can enlarge the stability region of E 0 .

Turing instability
Next, we will verify that system (2) can form steadily spatial patterns even if its firstderivative counterpart cannot develop any steady pattern. Firstly, we discuss the Turing instability of system (2). The unbalanced changes of phases, corresponding to Turing branches, are the transitions of system (2) from the uniform state to the oscillatory state [29,32]. After the process, the formed patterns are called Turing patterns. From (5) we obtain the necessary conditions for causing Turing instability when τ = 0: Conditions (8)-(10) indicate that system (2) is unstable for some perturbations to the wave number k. Thus, we obtain that det k (J) = 0 at the critical value. That is to say, Turing bifurcation occurs when Im(λ k ) = 0, Re(λ k ) = 0 at k = k c = 0 [19]. When Turing patterns come into being, the wave number k c satisfies Based on above discussions, we get the following result.
Lemma 3. If the parameters ϕ 1 and ϕ 2 are in the set Π 1 ∩ Π 2 and the following conditions hold then Turing instability in system (2) occurs when τ = 0.
When τ = 0, Turing instability occurs when λ → 0, namely, p 2 (k) + r → 0, and then we have the following Turing instability curve: Letting ϕ 2 (k) = ϕ 2 (k + 1) yields (11) and (12) hold . We set (ϕ 1 , ϕ 2 ) as free parameters, then project the parameter set Π 3 into a twodimensional space with respect to (ϕ 1 , ϕ 2 ). From Fig. 2 we observe that Π 1 ∩ Π 2 ∩ Π 3 corresponds with regions C (green) and D (yellow), where the Turing instability occurs. When the fraction order q = 1, the Turing instability region only consists of region C. That is to say, region D is the enlarged region, where Turing instability is caused by fractional derivatives and reaction-diffusion terms together.

Remark 1.
Usually, the curves of ϕ * 2 (ϕ 1 ) with k ∈ N may reduce the area of Turing instability region. In Fig. 2, regions C and D are exactly below the curve of ϕ * 2 (ϕ 1 ) with k = 2, and the Turing-Turing bifurcation point with k = 1 is on the left side of regions C and D. Thus, in this paper, the curves of ϕ * 2 (ϕ 1 ) with k ∈ N make no difference with the Turing instability region.

Numerical simulations
Numerical algorithms for reaction-diffusion systems are often complicated, which need tedious MATLAB programming. Garvie [6] proposed a semiimplicit (in time) finitedifference scheme to approximate the solutions of reaction-diffusion systems. The semiimplicit method means this algorithm involves approximations at the current time level t n and at the previous time level t n−1 . The corresponding algorithm is reduced to a sparse, banded and linear system of algebraic equations. In this section, we use Simulink to simplify the programming, which is more efficient and visually represents the interrelationships of neuron states. Simulink is graphical programming that is easy for nonprofessional scholars to program and realize the equivalent effect gained by MATLAB code. Besides, we can utilize user-defined MATLAB function blocks in Simulink, which takes full advantage of MATLAB code and Simulink.
Firstly, we give the following result about the Laplacian operator, which is used for the approximate calculation of reaction-diffusion terms. By means of the finite-difference scheme mentioned in [6,25], the Laplacian with respect to the concentration field w in the node (i, j) is calculated along the x and y directions simultaneously: where .
h x and h y in Eq. (13) are x and y grid spacings, respectively. Assume h = h x = h y and discretize Laplacian operation in a two-dimensional domain Ω: From (14) the two-dimensional Laplacian operator is Following the approach of the modified Adams-Bashforth-Moulton predictor-corrector scheme mentioned in [3], we derive the approximate numerical values of the fractional derivatives with delay. Denote ∆t as the time-step. Consider a uniform grid {t n = n∆t: n = −k, −k + 1, . . . , −1, 0, 1, . . . , N }, where k and N are integers such that k = τ /∆t and N = T /∆t. Discretize the initial condition Besides, the delayed terms can be rewritten as follows: where the discrete right-hand terms in system (2) are and the predictor terms arẽ In addition, a j,n+1 and b j,n+1 are defined as System (2) is simulated numerically in a 100 × 100 two-dimensional square region. Time step and space step are set as 0.005 and 0.5, respectively, which needs large amount of computation but ensures the accuracy of numerical simulations. Diffusion coefficients (d 1 , d 2 ) are chosen as (0.1, 1.6). Besides, c 1 = 2, c 2 = 4, a 1 = −4, a 2 = 2. b 1 and b 2 are chosen as free parameters. Simulink construction of system (2), as we can see in Fig. 3, consists of some blocks and lines. Each block can realize different functions and is connected with others by data flow lines. The main blocks are 2D Convolution and userdefined MATLAB function. The 2D convolution completes the approximate calculation of the reaction-diffusion terms by two-dimensional Laplacian operator L. The user-defined MATLAB function realizes the integration of the right-hand terms in system (2) by the above modified Adams-Bashforth-Moulton predictor-corrector scheme.
When fractional order q = 1 and τ = 0, system (2) becomes a first-derivative system. Activation function f (·) is selected as tanh(·), where b 1 = ϕ 1 and b 2 = ϕ 2 . Set (b 1 , b 2 ) as (4, 1) (marked by " * ") and (3.5, 3) (marked by "+"), respectively (as we can see in Fig. 2). Figures 4 and 5 are the two-dimension spatial states of neurons u(x, y, t) and v(x, y, t) at some moments, respectively. On account of diffusion, neuronic states not only vary in time, but also change in space. What is interesting is that spatial states take on particular shapes. Comparing the figures of u(x, y, t), v(x, y, t)( T = 50) with the figures of u(x, y, t), v(x, y, t) (T = 100) in Fig. 4, we observe that system (2) forms steady Turing patterns. Furthermore, when (b 1 , b 2 ) is (3.5, 3), we observe in Fig. 5 that as the time T increases from 50 to 100, the figures of u(x, y, t), v(x, y, t) (T = 50) are different from the figures of u(x, y, t), v(x, y, t) (T = 100), which means that under this condition, system (2) cannot induce any steady Turing pattern.
As we can observe in Fig. 6, system (2) with q = 0.79 can form steady spatial patterns. Comparing it with Fig. 5, we know that system (2) form steadily spatial patterns even if its first-derivative counterpart cannot develop any steady pattern. Hence, these patterns are induced by temporal fractional derivative and reaction-diffusion terms together, which corresponds with region D in Fig. 2. Furthermore, to investigate the relationship between fractional derivative and Turing pattern, we change the fractional order q as 0.72 or 0.87. From Figs. 7 and 8 we observe that the shape of Turing patterns varies with the fractional order q. Obviously, when q = 0.72, the patterns of u(x, y, t) and v(x, y, t) become various steady types. When q = 0.87, the stability performance of the patterns of u(x, y, t) and v(x, y, t) is a little weak. In Fig. 8, when the time T increases from 50 to 100, the center patterns change to some extent.
It rises a question that whether temporal fractional derivative has some relationship with the locally asymptotical stability region or Turing instability region with respect to (ϕ 1 , ϕ 2 ). From the left figure in Fig. 9 we obtain that the stability region of E 0 varies with fractional order q. As q increases from 0.72 to 0.79, the stability region of E 0 loses region D instead. Increasing q from 0.79 to 0.87 continuously, the stability region does not contain region C, namely, when q = 0.87, E 0 is stable when (ϕ 1 , ϕ 2 ) only stays in regions A and B. Besides, we wonder how the Turing instability region induced by fractional derivatives and reaction-diffusion terms together varies with fractional order q. From the right figure in Fig. 9, when q = 0.72, the Turing instability region induced by fractional derivatives and reaction-diffusion terms together consists of regions E, F and G.    The steadily spatial patterns of u(x, y, t) and v(x, y, t) for system (2) in two-dimension square domain with q = 0.72 and other parameters are similar to those in Fig. 6. Figure 8. The unsteadily spatial patterns of u(x, y, t) and v(x, y, t) for system (2) in two-dimension square domain with q = 0.87 and other parameters are similar to those in Fig. 6. As q increases, the corresponding Turing instability region shrinks to region E, that is to say, when (ϕ 1 , ϕ 2 ) stays in regions F and G, system (2) with q = 0.87 cannot form steady Turing patterns.
Remark 2. Due to that time delay widely exists in the neural network, we investigate the impact of transmission delay on Turing patterns. As we observe in Fig. 10, the patterns of u(x, y, t) and v(x, y, t) are more complicated than the patterns of those in Figs. 4-5 (a nondelay integer-order subsystem of system (2)) and Figs. 6-8 (a nondelay fractionalorder subsystem of system (2)). Compared with Fig. 6, as τ increase from 0 to 0.6, Turing patterns change unpredictably and become unstable.  Figure 10. The unsteadily spatial patterns of u(x, y, t) and v(x, y, t) for system (2) in two-dimension square domain with q = 0.79, τ = 0.6 and other parameters are similar to those in Fig. 6. Figure 11. The bifurcation diagrams of u(0, 0, ·) and v(0, 0, ·) with respect to the parameter τ ∈ [0, 1].
Remark 3. As shown in Fig. 11, the central points of u(x, y, t) and v(x, y, t) are selected to study the bifurcation diagrams of system (2) with respect to transmission delay τ . u(0, 0, ·) and v(0, 0, ·) stay in stable state when τ < τ 0 = 0.06 and oscillate greatly as τ exceeds τ 0 , where τ 0 can be seen as the critical value of Hopf bifurcation. Hence, as the transmission delay increases, Turing patterns lose its stability, and state variables of neurons vary in both time and space. Furthermore, when τ < τ 0 and fractional order is in a certain section, a stable Turing pattern may occur, which is induced by temporal fractional derivative and reaction-diffusion terms together; when τ > τ 0 , no matter how fractional order changes, the stable Turing pattern is hard to form.

Conclusion
In this paper, we proposed a delayed reaction-diffusion neural network with the Caputotype fractional derivative. The condition of Hopf bifurcation was obtained firstly through analyzing relevant characteristic equation, while the condition of Turing instability was derived following the pattern dynamics theory proposed by [23]. Especially, we found that fractional derivative has a certain relationship with Hopf bifurcation and Turing pattern.
To be specific, the fractional derivative could enlarge the stability region of homogeneous steady state E 0 concerning (ϕ 1 , ϕ 2 ), which means that when the first-derivative counterpart of system (2) is unstable at E 0 , but once introducing fractional derivative, E 0 may become stable again. Besides, as the fraction derivative decreases, it softens the Turing instability conditions of system (2). With the help of Simulink, we obtained the approximation solutions of system (2). Numerical simulations show that temporal fractional derivative contributes to Turing patterns in system (2), which accords with the previous theoretical analysis. Meanwhile, we found that the shape of the Turing pattern is related to the fraction order and the time delay. Pattern dynamics in neural networks with diffusion has been investigated in [4,14,29], but previous research does not consider fractional derivative. In our study, we illustrate the impact of fractional derivative on spatial and temporal dynamics for system (2) and obtain plenty of meaningful findings. [29] obtained the amplitude equations for system (1) and studied the selection of Turing patterns. Due to the existence of fractional derivative, it is difficult for us to obtain the amplitude equations of system (2). In our further research, we will try to improve some classic methods, such as multiple-scale analysis, to obtain the amplitude equations for fractional reaction-diffusion systems.
The spatial dynamics analysis of reaction-diffusion neural networks not only reveals some properties of integrated circuits that are used to build neural networks, where a Turing pattern corresponds to a steady state that leads to nonuniformly spatial oscillation, but also illustrates many biological phenomena such as normal neuron firing in the brain [18] or biological disorders such as fibrillation [24].