Edge of chaos in reaction diffusion CNN model

Abstract In this paper, we study the dynamics of a reaction-diffusion Cellular Nonlinear Network (RD-CNN) nodel in which the reaction term is represented by Brusselator cell. We investigate the RD-CNN dynamics by means of describing function method. Comparison with classical results for Brusselator equation is provided. Then we introduce a new RD-CNN model with memristor coupling, for which the edge of chaos regime in the parameter space is determined. Numerical simulations are presented for obtaining dynamic patterns in the RD-CNN model with memristor coupling.


Introduction
Nonlinear reaction-diffusion (RD) type of equations are widely used to describe phenomena in different fields, such as biology-Fisher model, Hodgkin-Huxley model and its simplification-FitzHugh-Nagumo nerve conduction model, physics -Sine-Gordon model, chemistry -Brusselator equation, et-cetera [1].
Semiconductor reaction-diffusion (RD) large scale circuits (LSI) implementing RD dynamics are considered widely in the literature [2,3]. These LSI, so called reaction-diffusion chips, are mostly designed by digital, analog, or mixed-signal complementary-metal-oxide-semiconductor (CMOS) circuits of Cellular Nonlinear Networks (CNN) and Cellular Automata (CA). In this paper we shall study a new reaction-diffusion CNN (RD CNN) model with memristors coupling. Cellular Nonlinear Networks (CNN) [2] are novel class of information processing systems constructed as an analog dynamic processor array which reflects just this property: the processing elements interact directly within a finite local neighborhood. Reaction-diffusion CNN [4] can virtually represent any feature of reaction-diffusion systems. It is known [2,4] that some autonomous CNN represent an excellent approximation to nonlinear partial differential equations (PDEs). The intrinsic space distributed topology makes the CNN able to produce real-time solutions of nonlinear PDEs. Generally, RD system can be described by: @u.x/ @t D g u r 2 u.
x/ C f v .u.x/; v.x//; where g u;v are the diffusion coefficients, f u;v .:/ is the reaction model.
In this paper, we shall consider three-molecule Brusselator [4] for the reaction term: where a; b are the reaction parameters. RD Brusselator model [5] is one of the first mathematical models used to explain the self-organization phenomena observed in chemical reactions of the reaction-diffusion type. Prigogine introduced dissipative structures which operate far from thermodynamic equilibrium. Two types of emergent behavior were found in a neighbourhood of a critical bifurcation boundary for the uncoupled cell: dynamic patterns or stationary patterns. Here, we shall predict a sharper precise domain in the cell parameter space where emergent behaviour appears, called edge of chaos.
In Section 2, we introduce the RD CNN model with Brusselator cell (RD-BCNN) in the reaction term. Section 3 deals with its dynamics by applying describing function method. Comparison with the classical results [5] will be presented as well. In Section 4 we shall study the dynamics of RD-BCNN model with memristors coupling. We shall find the domain of cell parameters in which RD-BCNN model will operate far from thermodynamic equilibrium.

RD CNN model with Brusselator cell
Cellular Nonlinear Networks (CNN) are analogue dynamic processor arrays, which are made of cells. Let us consider a two-dimensional grid with 3 3 neighborhood system shown on Figure 1. One of the key features of CNN is that the individual cells are nonlinear dynamical systems, but the coupling between them is linear. Roughly speaking, one could say that these arrays are nonlinear but have a linear spatial structure, which makes the use of techniques for their investigation common in engineering or physics attractive.
We will give general definition of a CNN which follows the original one [2]: 1. An M M cellular nonlinear network is defined mathematically by four specifications: 1) CNN cell dynamics; 2) CNN synaptic law which represents the interactions (spatial coupling) within the neighbor cells; 3) Boundary conditions; 4) Initial conditions.
In terms of the definition we can present the dynamical systems describing CNN. For general CNN whose cells are made of time-invariant circuit elements, each cell C.ij / is characterized by its CNN cell dynamics : where x ij 2 R m and u ij is usually a scalar. In most cases, the interactions (spatial coupling) with the neighbor cell C.i C k; j C l/ are specified by a CNN synaptic law: The first term A ij;kl x i Ck;j Cl of (4) is simply a linear feedback of the states of the neighborhood nodes. The second term provides an arbitrary nonlinear coupling, and the third term accounts for the contributions from the external inputs of each neighbor cell that is located in the N r neighborhood. CNN is an excellent candidate for both analog and digital applications because of its structural simplicity, relative ease of fabrication, inherent speed and design flexibility. Many methods used in image processing and pattern recognition can be easily implemented by CNN approach.
There are several ways to approximate the Laplacian operator r 2 in discrete space by a CNN synaptic law with an appropriate A-template [4]. We choose the following template for the two-dimensional reaction-diffusion CNN model: Then our CNN model for the RD Brusselator equations (1) and (2) with A-template can be written in the following form (RD-BCNN): The RD Brusselator equations (1) and (2) is well known in chemical kinetics as an ideal system for studying dissipative structures. In some sense, this system behaves like harmonic oscillator. In [5] it is shown that the dynamics of RD Brusselator equation has global periodicity which depends on the initial and boundary conditions. We shall sdudy the dynamical behavior of the RD-BCNN model (6) in the next section.

Dynamics of RD-BCNN
For the analysis of this CNN model (6) we shall apply the describing function method. In this case the state variable u kl .t/ are functions depending on three arguments -two discrete space arguments k; l, and one continuous time t . We use the following double Fourier transform F .s; z 1 ; z 2 / of functions f kl .t /, continuous in time and discrete in space: Applying the above transform to our CNN model (6) we obtain: where s D i! 0 , z 1 D exp.i 1 /, z 2 D exp.i 2 /, ! 0 , 1 , 2 are temporal and two spatial frequencies, respectively and i D p 1 is the identity. According to the describing function method [6], the transfer function H.s; z 1 ; z 2 / can be obtained from the above equation as: Since we are looking for eventual periodic solutions, we shall suppose that the state has the following form: where the temporal frequency is ! 0 D 2 T 0 , T 0 > 0 is the minimal period. If we take the periodic boundary conditions for our CNN model (6), making the array circular, we get: According to the describing function method [6], we shall approximate the output v kl .t / by the fundamental component of its Fourier expansion: with where f is the piecewise-linear sigmoid output function f .x/ D 1 2 .jx C 1j jx 1j/. So we obtain: V m 0 D 2U m 0 . We shall express now the transfer function H from (9) in terms of ! 0 , 1 , 2 : According to (9) and (14) the following constraints on its real and imaginary parts of the transfer function hold: If, for a given value of 1 , 0 Ä 1 < 2 , we can find a solution .! 0 ; 2 ; U m 0 ; V m 0 / to these two equations, the describing function method predicts the existence of a periodic solution with an amplitude approximately U m 0 and a minimal period T 0 D 2 ! 0 [6]. Therefore from (11), (15) and (16) we get the following system: From the second equation of (17) we obtain B D 1 D 1 D 2 . Then combining with the first equation of (17) we get: Solving the third equation of (17) with respect to ! 0 we obtain: Therefore the following theorem has been proved: We shall make a comparison with the classical results for RD Brusselator equations (1) and (2). Let us linearized the RD Brusselator equation (1), (2) using the following change of variables: Then the linearized system with respect to the variables Q u, Q v can be written as: where the operator L is parabolic of the following form: Then the solutions of the linearized RD Brusselator equation are: si n m Ä n ; m D 1; 2; : : : ; 0 Ä Ä Ä n; n is the dimension of the linearized system. According to [5] for Dirichlet boundary conditions the RD Brusselator equation will have periodic solutions if the diffusion coefficients satisfy the following inequality: Moreover, if D 1 D D 2 D D, we obtain D < an 2 2 . In other words, if the values of the diffusion coefficients are close to each other then periodic solutions arise for the RD Brusselator equations (1) and (2).
There are many interesting applications in biology, in which chemical reactions can be produced in the conditions of closed geometry. For example, a surface of a sphere, which is the cell membrane. Then we have periodic boundary conditions of the type: In other words nonhomogeneous solutions with a space periodicity arise. Now if we compare the obtained results for RD-BCNN model (6) in this section with the classical ones (20), we can see the similarity of the periodic solutions.

RD-BCNN with memristor coupling
Let us consider equations (1) and (2). When u.x/ and v.x/ are represented by voltages on the RD hardware, the gradient (diffusion terms in the RD model) is represented by linear resistors. We discretize first equation of (1) in the following way: where s is the spatial index, x is the discrete step in space, terms g u .u s 1 u s / and g u .u sC1 u s / represent respectively current flowing into the sth node from .s 1/th and .s C1/th nodes via two resistors whose conductance is represented by g u . The spatial Laplacian r 2 in equation (1) can be approximated as Here, we introduce the memristor model [7], in which the resistors are replaced with memristors. The resulting point dynamics is given by where g u .:/ denotes the monotonically increasing function defined by: where the right-hand side represents the current of the memristors [7], Á 1 denotes the polarity coefficient - In this section we shall identify the values of the cell parameters for which the interconnected RD system may exhibit complexity. The necessary condition for a nonconservative system to exhibit complexity is to have its cell locally active [3]. The theory which will be presented below offers a constructive analytical method for uncovering local activity. The precisely defined parameter domain in which the model can exhibit complex behavior is called the edge of chaos. The RD system operating near the edge of chaos, where the cells are not locally active, is also linearly asymptotically stable. In particular, constructive and explicit mathematical inequalities can be obtained for identifying the region in the CNN parameter space where complexity phenomena may emerge, as well as for localizing in further into a relatively small parameter domain called edge of chaos where the potential for emergency is maximized. By restricting the cell parameter space to the local activity domain, a major reduction in the computing time required by the parameter search algorithms is achieved [3].
We shall apply the theory of local activity in order to study dynamic behaviour of RD-CNN model with memristors coupling in the following form: where D p .w L;R / denotes the monotonically increasing function defined as (26) Then the memristive dynamics is defined as in (24).
We develop the following constructive algorithm for determining the edge of chaos domain.
1. Map the RD-BCNN with memristor coupling into its discrete-space version: 2. Find the equilibrium points of RD-BCNN model (27). According to the theory of dynamical systems equilibrium points u , v are these for which In general, system (28) may have one, two or three real roots and these roots are functions of the cell parameters a; b. In the case the Brusselator cell we have a unique equilibrium point E 1 D .
3. Calculate the cell coefficients of the Jacobian matrix of (28) about the system equilibrium point E 1 , namely, 4. Calculate the trace T r.E 1 / and the determinant .E 1 / of the Jacobian matrix for the unique equilibrium point.
Remark. The importance of the circuit model is not only in the fact that we have a convenient physical implementation, but also in the fact that well-known results from classic circuit theory can be used to justify the cells' local activity. In this sense, if there is at least one equilibrium point for which the circuit model of the cell acts like a source of "small signal", i.e. if the cell is capable of injecting a net small-signal average power into the passive resistive grids, then the cell is said to be locally active.

Edge of chaos
In the literature, the so-called edge of chaos (EC) means a region in the parameter space of a dynamical system where complex phenomena and information processing can emerge. We shall try to define more precisely this phenomena till now known only via empirical examples.
Definition 4.2. RD-BCNN with memristor coupling is said to be operating on the edge of chaos EC iff there is at least one equilibrium point which is both locally active and stable.
In the case of Brusselator cell we prove the following theorem: Proof. The Jacobian matrix A of the system (28) is: a 11 a 12 a 21 a 22 where a 11 D @f 1 @u kl D b 1, a 12 D @f 1 @v kl D a 2 , a 21 D @f 2 @u kl D b, a 22 D @f 2 @v kl D a 2 .
Since we found the unique equilibrium point of the RD-BCNN to be E 1 D .a; b a /, then .
After simulations, we obtain the following edge of chaos region for the RD-BCNN with memristor coupling: The local activity domain and the edge of chaos domain are shown on Fig. 2 for the following parameter set: a 2 OE 5; 5 and b 2 OE 0:4; 0:4. Most of the cell parameters points, which were found capable of self-organization are located nearby the bifurcation boundary separating the stable region from the unstable region. This simulation provides very useful information concerning the potential of the cell to exhibit complexity. All simulations obtained in the paper use MATCNN software [8]. Now, we shall present a short discussion of the dynamic nonhomogeneous patterns which have emerged from the RD-BCNN with memristor coupling for several cell parameter points chosen arbitrarily within, or nearby, the edge of chaos domain. After extensive computer simulations we obtain the following results: The dynamic patterns presented in Figure 3 are associated with cell parameter points located within the EC region, which corresponds to the locally active and stable domain. Dynamic behaviour in this case will correspond to cells oscillating independently of each other with uncorrelated phases. In a visual representation of the type shown in Figs. 3 such a nonsynergetic behavior would be seen as a change in the shape of the patterns during the temporal evolution. At a glance, one can see that this is clearly in all of these cases we have obtained an emergent, complex behavior.
During the simulation procedure the initial state of all Brusselator cells is set to be inactive state. After stimulating the center node, the excitable waves propagate outwards, resulting in the aggregation of patterns. Figure 3 shows the surface wave patterns of RD-BCNN with memristor coupling. Although the center was stimulated, wave propagates asymmetrically due to memristive effects with same polarity over the medium. The velocity of waves propagation was differed depending on the direction of the wave propagation.

Conclusions
In this paper, we studied RD CNN model with Brusselator cell. We started by studying the dynamics of RD-BCNN, applying describing function method and we compare the obtained results with the classical ones for RD Brusselator equation [5]. Then we introduced memristor dynamics in our RD-BCNN and determined the edge of chaos domain of the parameter set in which the model exhibits complex dynamical behavior. We also provided extensive computer simulations of the dynamic patterns which can be obtained in the edge of chaos domain.
Among the demonstrated behaviour (Fig. 3), the nonuniform spatial pattern generation is applied to investigation on detecting global motion of excitable waves, because the memristors [7] accumulated the directions of waves on the media, which resulted in detecting majority of wave directions at every spatial point. This application is not limited to the analysis of RD systems, but the idea can be transferred to analyzing much more complex systems like brain networks, social networks and so on.
Recent laboratory observations suggesting that chaotic regimes may in fact represent the ground state of central nervous system point to the intriguing possibility of exploiting and controlling chaos for future scientific and engineering applications [3].