Edge of chaos in reaction diffusion CNN model

Angela Slavova; Ronald Tetzlaff

doi:10.1515/math-2017-0002

Article Open Access

Edge of chaos in reaction diffusion CNN model

Angela Slavova and Ronald Tetzlaff

Published/Copyright: February 2, 2017

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

$Open Mathematics$

From the journal Open Mathematics Volume 15 Issue 1

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.

Keywords: Cellular Nonlinear Networks; Reaction-diffusion Brusselator equation; Memristor; Edge of chaos

MSC 2010: 92B20; 35Q80; 65P20

1 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:

(1) ∂u(x)∂t=gu∇2u(x)+fu(u(x),v(x));∂v(x)∂t=gv∇2v(x)+fv(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:

(2) fu(u(x),v(x))=a−(b+1)u+vu2,fv(u(x),v(x))=bu−vu2,

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.

2 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 x 3 neighborhood system shown on Figure 1.

$Fig 1 3 x 3 neighborhood CNN$

Fig 1

3 x 3 neighborhood CNN

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]:

Definition 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 :

(3) x ˙ i j = − g ( x i j , u i j , I i j s ) ,

where x_ij ∊R^m and u_ij is usually a scalar. In most cases, the interactions (spatial coupling) with the neighbor cell C(i + k, j + l)are specified by a CNN synaptic law:

(4) I i j s = A i j , k l X i + k , j + l + A ~ i j , k l ∗ f k l ( x i j ; x i + k , j + l ) + B ~ i j , k l ∗ u i + k , j + l ( t )

The first term A_ij,klx_{i+k,j +l} 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 ∇² 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:

(5) A = 0 1 0 1 − 4 1 0 1 0 .

Then our CNN model for the RD Brusselator equations (1) and (2) with A-template can be written in the following form (RD-BCNN):

(6) d u k l d t = a − ( b + 1 ) u k l + u k l 2 v k l + D 1 [ u k + 1 l + u k − 1 l + u k l + 1 + u k l − 1 − 4 u k l ] , d v k l d t = b u k l − u k l 2 v k l + D 2 [ v k + 1 l + v k − 1 l + v k l + 1 + v k l − 1 − 4 v k l ] , 1 ≤ k ≤ M , 1 ≤ l ≤ M .

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.

3 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₁, z₂)of functions ƒ_kl (t),continuous in time and discrete in space:

(7) F ( s , z 1 , z 2 ) = ∑ k = − ∞ ∞ z 1 − k ∑ k = − ∞ ∞ z 2 − l ∫ − ∞ ∞ f k l ( t ) e x p ( − s t ) d t .

Applying the above transform to our CNN model (6) we obtain:

(8) s U = − ( b + 1 ) U + U 2 V + D 1 [ z 1 U + z 1 − 1 U + z 2 U + z 2 − 1 U − 4 U ] , s V = b U − U 2 V + D 2 [ z 1 V + z 1 − 1 V + z 2 V + z 2 − 1 V − 4 V ] ,

where s = iω₀, z₁ = exp(iΩ₁), z₂ = exp(iΩ₂), ω₀, Ω₁, Ω₂ are temporal and two spatial frequencies, respectively and i = − 1 is the identity.

According to the describing function method [6], the transfer function H(s, z₁, z₂) can be obtained from the above equation as:

(9) H ( s , z 1 , z 2 ) = D 2 ( z 1 + z 12 − 1 + z + z 2 − 1 − 4 ) − s s + 1 − D 1 ( z + z 12 − 1 + z + z 2 − 1 − 4 ) = U V .

Since we are looking for eventual periodic solutions, we shall suppose that the state has the following form:

(10) u k l ( t ) = U m 0 s i n ( ω 0 t + k Ω 1 + l Ω 2 ) ,

where the temporal frequency is ω 0 = 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:

(11) Ω 1 + Ω 2 = 2 K π n , 0 ≤ K ≤ n − 1 , n = M . M

According to the describing function method [6], we shall approximate the output υkl(t) by the fundamental component of its Fourier expansion:

(12) v k l ( t ) ≈ V m 0 s i n ( ω 0 t + k Ω 1 + l Ω 2 ) ,

with

(13) V m 0 = 1 π ∫ − π π f ( U m 0 s i n ψ ) d ψ ,

where ƒ is the piecewise-linear sigmoid output function f ( x ) = 1 2 ( | x + 1 | − | x − 1 | ) . So we obtain: V_mo = 2U_mo. We shall express now the transfer function Η from (9) in terms of ω₀, Ω ₁, Ω₂:

(14) H Ω 1 , Ω 2 ( ω 0 ) = D 2 ( 2 c o s Ω 1 + 2 c o s Ω 2 − 4 ) − i ω 0 i ω 0 + 1 − D 1 ( 2 c o s Ω 1 + 2 c o s Ω 2 − 4 ) .

According to (9) and (14) the following constraints on its real and imaginary parts of the transfer function hold:

(15) R ( H ) = D 2 B − D 2 D 1 B 2 − ω 0 2 ( 1 − D 1 B ) 2 + ω 0 2 = U m 0 V m 0 ;

(16) J ( H ) = ( D 1 B − 1 − D 2 B ) ω 0 ( 1 − D 1 B ) 2 + ω 0 2 = 0 ,

where Β = 2cosΩ₁ + 2cosΩ₂ - 4.

If, for a given value of Ω1, 0 ≤ Ω1 < 2π, we can find a solution (ω₀, Ω₂, U_mo, V_mo ) to these two equations, the describing function method predicts the existence of a periodic solution with an amplitude approximately U_mo and a minimal period T 0 = 2 π ω 0 [6]. Therefore from (11), (15) and (16) we get the following system:

(17) Ω1+Ω2=2Kπn;D1B−1−D2B=0;D2B−D2D1B2−ω02(1−D1B)2+ω02=12.

From the second equation of (17) we obtain B = 1 D 1 − D 2 . Then combining with the first equation of (17) we get:

Ω1=Kπn+arccos[4(D1−D2)+14cosKπn(D1−D2)],Ω2=Kπn−arccos[4(D1−D2)+14cosKπn(D1−D2)],0≤K≤n−1.

Solving the third equation of (17) with respect to ω₀ we obtain:

ω 0 = B ( D 2 + 4 D 1 ) − B 2 ( D 2 D 1 + 2 D 1 2 ) − 2 3 .

Therefore the following theorem has been proved:

Theorem 3.1

CNN model (6) of the RD Brusselator equations (1) and (2) with circular array of η = M.M cells has periodic solution with period T 0 = 2 π ω 0 and amplitude U m 0 for all Ω 1 + Ω 2 = 2 K π n , 0 ≤ K ≤ n − 1.

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:

u=a+u~,v=ba+v~.

Then the linearized system with respect to the variables u͂, v͂ can be written as:

∂ ∂ t u ~ v ~ = L u ~ v ~ ,

where the operator L is parabolic of the following form:

L = b − 1 + D 1 ∇ 2 a 2 − b − a 2 + D 2 ∇ 2 .

Then the solutions of the linearized RD Brusselator equation are:

(18) u ~ v ~ = ∑ m d m e ω m t c 1 c 2 s i n m π κ n , m = 1 , 2 , … , 0 ≤ κ ≤ 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:

(19) D 2 − D 1 < n 2 π 2 ,

Moreover, if D₁ = D₂ = D, we obtain D < a n 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:

u ~ ( 0 ) = u ~ ( n ) ; v ~ ( 0 ) = v ~ ( n ) .

In this case the periodic solutions of the RD Brusselator equation are of the following form:

(20) u ~ v ~ = c 1 c 2 e ω m t 2 π i q κ n e , q = 0 , 1 , …

In other words nonhomogeneous solutions with a space periodicity arise. Now if we compare the obtained results for RD-BCNN model (6)s in this section with the classical ones (20), we can see the similarity of the periodic solutions.

4 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:

(21) d u s d t = g u ( u s − 1 − u s ) + g u ( u s + 1 − u s ) Δ x 2 + f u ( . ) ,

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_s+1- u_s)represent respectively current flowing into the sth node from (s – l)th and (s + l)th nodes via two resistors whose conductance is represented by g_u. The spatial Laplacian V² in equation (1) can be approximated as

∇ 2 u ( x ) = u s − 1 + u s + 1 − 2 u s ∇ x 2 .

Here, we introduce the memristor model [7], in which the resistors are replaced with memristors. The resulting point dynamics is given by

(22) d u s d t = g u ( w s L ) ( u s − 1 − u s ) + g u ( w s R ) ( u s + 1 − u s ) Δ x 2 + f u . ; d v j d t = f v . ,

where g_u(.)denotes the monotonically increasing function defined by:

(23) g u w s L , R = g m i n + g m a x − g m i n 1 1 + e − β w s L , R ,

where β denotes the gain, g_min and g_max denote the minimum and maximum coupling strengths, respectively, w S L . R denote the variables for determining the coupling strength ( L - leftward, R - rightward). Finally, we introduce the following memristive dynamics for w S L , R :

(24) τ d w s L , R d t = g u ( w s L , R ) . η 1 . ( u s − 1 − u s ) ,

where the right-hand side represents the current of the memristors [7], n₁denotes the polarity coefficient -n₁ η 1 = + 1 : w s L , η 1 = − 1 : w s R .

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:

(25) dukldt=f1(ukl,vkl)+D1(wL,R)(uk−1l+uk+1l+ukl−1+ukl+1−4ukl);dvkldt=f2(ukl,vkl)+D2(wL,R)(vk−1l+vk+1l+vkl−1+vkl+1−4vkl),

where D_p(w^L,R)denotes the monotonically increasing function defined as

(26) D p ( w L , R ) = D m i n ( D m a x − D m i n ) . 1 1 − e − β w p L , R , p = 1 , 2 ,

f₁(u_kl,v_kl) = a-(b + 1)u_kl+v_klu²_kl , f₂(u_kl,v_kl) = bu_kl-v_klu²_kl. 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:

(27) d u k l d t = f 1 ( u k l , v k l ) + D 1 ( w k l L , R ) ( u k − 1 l + u k + 1 l + u k l − 1 + u k l + 1 − 4 u k l ) ; d v k l d t = f 2 ( u k l , v k l ) + D 2 ( w k l L , R ) ( v k − 1 l + v k + 1 l + v k l − 1 + v k l + 1 − 4 v k l ) ,

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

(28) f 1 ( u ∗ , v ∗ ) = 0 ; f 2 ( u ∗ , v ∗ ) = 0.

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 = ( u 1 ∗ , v 1 ∗ ) = ( a , b a ) .

3. Calculate the cell coefficients of the Jacobian matrix of (28) about the system equilibrium point E₁, namely,

a 11 ( E 1 ) = b − 1 ; a 12 ( E 1 ) = a 2 ; a 21 ( E 1 ) = − b ; a 22 ( E 1 ) = − a 2 .

4. Calculate the trace Tr(E₁)and the determinant Δ(E₁)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.

Definition 4.1

Stable and locally active region SLAR(E₁) at the equilibrium point E₁ for the RD-BCNN with memristor coupling model is such that

a 22 > 0 o r 4 a 11 a 22 < ( a 12 + a 21 ) 2

and

T r < 0 a n d Δ > 0.

5. 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:

Theorem 4.3

RD-BCNN with memristor coupling (27) is operating in the EC regime iff the following conditions for the parameters are satisfied: –4(b– 1)a² < (a² – b)² and b < a² + 1. In this parameter set the unique equilibrium point Ε₁ is both locally active and stable.

Proof. The Jacobian matrix A of the system (28) is:

(29) A = a 11 a 12 a 21 a 22 ,

where a 11 = ∂ f 1 ∂ u k l = b − 1 , a 12 = ∂ f 1 ∂ v k l = a 2 , a 21 = ∂ f 2 ∂ u k l = − b , a 22 = ∂ f 2 ∂ v k l = − a 2 .

Since we found the unique equilibrium point of the RD-BCNN to be E 1 = ( a , b a ) , then

(30) T r ( E 1 ) = − a 2 + b − 1 ,

(31) Δ ( E 1 ) = a 2 .

Now we shall calculate SLAR(E₁)according to Definition 4.1: SLAR(E₁): -4(b-1)a² < (a²-b²)²and b < a²+ 1. Therefore, according to Definition 4.2, the parameter set in which the RD-BCNN with memristor coupling is operating in edge of chaos regime is: –4(b – 1)a² < (a² – b²)²and b < a²+ 1. Theorem is proved.

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 ∊ [–5, 5] and b ∊ [–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].

$Fig 2 Edge of chaos region in RD-BCNN with memristor coupling$

Fig 2

Edge of chaos region in RD-BCNN with memristor coupling

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.

$Fig 3 Dynamic patterns in RD-BCNN with memristor coupling$

Fig 3

Dynamic patterns in RD-BCNN with memristor coupling

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.

5 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].

Acknowledgement

This first author is supported by Fellowship at Technical University of Dresden.

References

[1] Britton N.F., Reaction-Diffusion Equations and Their Applications to Biology, New York: Academic, 1986Search in Google Scholar

[2] Chua L.O., Yang L., Cellular neural networks: Theory, IEEE Trans. CAS, 1988, 35, 1257-127210.1109/ISCAS.1988.15089Search in Google Scholar

[3] Chua L.O., Local activity is the origin of complexity, International Journal of Bifurcation and Chaos, 2005, 15:11, 3435-345610.1109/CNNA.1998.685319Search in Google Scholar

[4] Slavova A., Cellular Neural Networks:Dynamics and Modeling, Kluwer Academic Publishers, 200310.1007/978-94-017-0261-4Search in Google Scholar

[5] Prigogine I., From Being to Becoming (Freeman, San Francisco), 1980Search in Google Scholar

[6] Mees A.I., Dynamics of feedback systems, London, Wiley, 1981Search in Google Scholar

[7] Chua L.O., Memristor: the missing circuit element, IEEE Trans. On Circuit Theory, 1971, 18:5, 507-51910.1109/TCT.1971.1083337Search in Google Scholar

[8] CNN Software Library, Ver1 .1, Analogical and Neural Computing Laboratory, Budapest, 2000Search in Google Scholar

Received: 2016-10-10

Accepted: 2016-12-20

Published Online: 2017-2-2

This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 3.0 License.