In this chapter we present mathematical study of memristor systems. More precisely, we apply local activity theory in order to determine the edge of chaos regime in reaction-diffusion memristor cellular nanoscale networks (RD-MCNN) and in memristor hysteresis CNN (M-HCNN). First we give an overview of mathematical models of memristors, CNN and complexity. Then we consider the above mentioned two models and we develop constructive algorithm for determination of edge of chaos in them. Based on these algorithms numerical simulations are provided. Two applications of M-HCNN model in image processing are presented.
- cellular nanoscale networks
- reaction-diffusion systems
- edge of chaos
- image processing
Memristors form an important emerging technology for memory and neuromorphic computing applications (Figure 1). Chua has developed the fundamentals of the memistor framework nearly 40 years ago . Since then, the industry has been engaged in the search for novel materials and technologies of these nano-structures .
Mathematical models of the complex dynamics which can be exhibited by nano-devices is presented in . General and simple models are very important in the investigations of nonlinear dynamics in memristors . Such models of memristor-based circuits are presented in [5, 6]. In order to develop novel hybrid  hardware architectures combining memory storage and data processing in the same physical location and at the same time  to explain the behavior of biological systems  new accurate mathematical models need to be introduced (Figure 2). Although several physical models [10, 11, 12, 13] have been derived in order to study phenomena characterizing these nano-devices, a circuit theoretic-based mathematical treatment allowing the development of memristor circuits is still restricted to few cases. Most of these mathematical models have been studied in . The merit for the first model of a titanium dioxide-based nano-structure may be ascribed to Williams . This model is simple not specifying boundary conditions in the state equation. Literature was later enriched with a number of more complex models taking into account nonlinear effects on ionic transport and defining behavior at boundaries. The model proposed by Joglekar and Wolf  may allow for single-valued state-flux characteristics only, under any sign-varying periodic input with zero mean. By contrast, only multi-valued state-flux characteristics may be reproduced by Biolek’s model  under the same type of excitation. A comparison of the dynamic behavior of titanium dioxide memristor circuits assuming different boundary models including the BCM model is given in . The results underline the sensitivity of these nonlinear circuits to the modeling accuracy. In the following, models have been proposed also for memristors based on other materials, e.g., for a tantalum oxide element [13, 19] and for a niobium dioxide memristors . Nevertheless, recent investigations  uncover numerical problems occurring in the numerical solution of the strongly nonlinear differential-algebraic equations or in a SPICE simulation  of these memristors.
Computing with memory is one of the main properties of neural networks, in which the performance is within localized memory storage. Therefore, they can provide high density analogue storage and can be integrated locally with computing elements which is the main advantage of the network learning algorithms. The memristive array works as a conventional memory, i.e., the weigh values can be calculated outside the array and can be programmed to the correct addresses. The algorithms are flexible because they can be implemented with an appropriate design. Local and non-local learning algorithms can be implemented in a straightforward way. The disadvantage is that they cannot have high connectivity and internal dynamics.
Cellular Nonlinear/Nanoscale Networks (CNN) have been introduced in 1988 by Chua and Yang [20, 21] as a new class of information processing systems which show important potential applications (Figure 3). By endowing a single CNN cell with local analog and logic memory, some communication circuitry and further units, the CNN Universal Machine (CNN-UM) has been invented by Roska [22, 23, 24]. Analog CNN-UM chip hardware implementations have been developed some time ago. A CNN-UM chip represents a parallel computer with stored programmability allowing real-time processing of multivariate data. CNN have very promising applications in image processing and pattern recognition [22, 25, 26, 27]. Although, recently realized systems, e.g., the EyeRis 1.3 system, the MIPA4k, and SCAMP-5, are characterized as sensor-processor systems for high speed vision by reaching frames rates more than 20 kHz, their low resolution (e.g., 176 × 144 pixel in the EyeRis 1.3 system) limits the applicability to certain problems in practice only. Since the cell size cannot be decreased considerably in conventional CMOS technology, nano-elements will play an important role in future CNN-UM chip realizations. Especially, memristors  which are considered for synaptic connections in first realizations , will play an important role for the realization of future CNN-UM sensor-processor systems by taking their rich dynamical behavior into account. However, a deep mathematical treatment of CNN with memristors, briefly called memristor CNN in the following, has not been provided so far. Especially, the derivation of methods allowing the determination of the parameter space of a memristor CNN showing emergent complex behavior, is being essentially important in the development of CNN-based computational methods.
Many scientists have struggled to uncover the elusive origin of “complexity” and its many equivalent jargons, such as emergence, self-organization, synergetics, collective behaviors, non-equilibrium phenomena, etc. [22, 30, 31, 32]. In his works, Schrödinger  defines the exchange of energy as a necessary condition for complexity of open systems. Prigogine  states the instability of the homogeneous systems as a new principle of nature, whereas Turing finds the origin of morphogenesis in symmetry breaking. In  Smale considers a reaction-diffusion system which has properties such that it makes the Turing interacting system to oscillate. 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.
Among other things, a mathematical proof is given in [22, 30, 33, 34] showing that none of the complexity related jargons cited above can explain emergent complex behavior in reaction-diffusion system without introducing local activity. The theory of local activity offers a constructive method for determination of complexity. We shall propose algorithm for hysteresis CNN which defines the domain of the cell parameters in which the system is capable of exhibiting complexity. The main advantage of local activity theory is that the complex behavior of reaction-diffusion system can be explained in a rigorous way by explicit mathematical formulas determining a small subset of locally active parameters’ region called edge of chaos. Cell kinetic equations which are locally active can exhibit limit cycles or chaos if the cells are uncoupled from each other by letting all diffusion coefficients to be zero. In this case complex spatiotemporal phenomena arise, such as spatiotemporal chaos or scroll waves.
In particular, constructive and explicit mathematical inequalities can be obtained for identifying that region in the CNN parameter space. 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 [33, 34].
In the next sections we shall present two mathematical models of memristor CNN. First one is reaction-diffusion memristor CNN, and the second one is memristor hysteresis CNN. We shall derive algorithm for determination of edge of chaos regime in these models based on local activity theory .
2. Reaction-diffusion memristor CNN (RD-MCNN)
Nonlinear reaction-diffusion types of equations are widely used to describe phenomena in different fields. We shall determine for reaction-diffusion CNN the domain of the cell parameters in which the cells are locally active and therefore they can exhibit complex behavior. Edge of chaos (EC) is associated with a region of parameter space in which complex phenomena and thus information processing can appear.
In this section the principle of local activity will be applied in studying complex behavior of reaction-diffusion CNN with memristor synapses (RD-MCNN). Semiconductor reaction-diffusion (RD) large scale circuits (LSI) implementing RD dynamics, called reaction-diffusion chips, are mostly designed by digital, analog, or mixed-signal complementary-metal-oxide-semiconductor (CMOS) circuits of CNN and cellular automata (CA).
In our model each cell will be arranged on a two-dimensional square grid and will be connected to adjacent cells through coupling devices that mimic 2-D spatial diffusion and transmit the cell’s state to its neighboring cells, as in conventional CNN.
We shall consider a discrete medium of identical cells which interact locally and therefore the homogeneous medium exhibits a non-homogeneous static or spatiotemporal patterns under homogeneous initial and boundary conditions. The theory of local activity will be formulated mathematically and implemented in circuit models. We shall start with reaction-diffusion CNN equations as a special class of spatially extended dynamical systems and we shall define the principle of local activity which will not be based on observations but on rigorous mathematical analysis.
2.1. Definition of local activity for reaction-diffusion equations
We call reaction-diffusion CNN equations locally active if and only if their cells are locally active at some cell equilibrium points .
In  principle of local activity is explained with the assumption of zero energy at the zero time. Therefore we can say that the cell is acting as a source of small signal energy and it is able to give rise to an initially very small input signal to a larger energy signal.
From mathematical point of view the signal should be very small in order to take the linear terms of Taylor series expansion of the cell model. In this way we can derive explicit analytical inequalities for the cell to be locally active at some equilibrium points in which the Taylor series expansion is computed. In other words, we can say that complex behavior of cells arises from infinitesimal small perturbations.
where are the diffusion coefficients, state for the reaction model. In (1)u(x) and v(x) are represented by voltages on the RD hardware, and the gradient is represented by linear resistors.
Let us discretize first equation of (1) in space:
where j is the spatial index, ∆x is the discrete step in space, terms and represent respectively current flowing into the jth node from (j − 1)th and (j + 1)th nodes via two resistors whose conductance is represented by .
We consider the memristor model , in which the resistors are replaced with memristors:
where the voltage across the memristor is v, the current of the memristor is i, the nominal internal state of the memristor corresponding to the charge flow is w, and the monotonically non decreasing function is when w is increasing.
We shall replace resistors for diffusion in analog RD LSIs with memristors (Figure 4).
Then we obtain the resulting point dynamics
where is the monotonically increasing function defined by:
where is the gain, and are the minimum and maximum coupling strengths, respectively, and denotes the variables for determining the coupling strength (l—leftward, r—rightward).
We introduce the following memristive dynamics for :
where the right-hand side represents the current of the memristors in (2), denotes the polarity coefficient—= + 1: , = −1: .
In the next we shall study RD-MCNN model of FitzHugh-Nagumo system.
Simplification of Hodgkin-Huxley model (Figure 5) can be given by FitzHugh-Nagumo system consisting of two coupled partial differential equations with two diffusion coefficients. In generally it describes the electrical potential across cell membrane when the flow of ionic channels is changed. It also can be presented as the model of electrical waves of the heart.
In this section we shall present the following FitzHugh-Nagumo system with two diffusion terms:
f(u) is monotonically non decreasing function
The parameters ε, a, b are physical parameters, related to pressures of oxide, carbon oxide and to the temperature.
We discretize spatially system (7) and the resulting point dynamics are given as:
where denotes the monotonically increasing function defined as
Then the memristive dynamics is defined as in (6).
We map memristive FitzHugh-Nagumo system into the associated FitzHugh-Nagumo CNN model:
We find the equilibrium points of (11). According to the theory of dynamical systems equilibrium points are these for which:
System (12) may have one, two or three real roots . In general, these roots are functions of the cell parameters a, b, ε.
We calculate the four cell coefficients of the Jacobian matrix of (12) about each system equilibrium point .
Then we calculate the trace ) and the determinant of the Jacobian matrix of (12) for each equilibrium point.
We define locally active region for each equilibrium point :
Additional stability condition in this case is:
It can be shown that this is the only region which corresponds to locally asymptotically stable equilibrium points of our model.
We define the stable and locally active region
In our particular case, we have three equilibrium points .
Then we check the conditions for the local activity and stability of the equilibrium points. The result is that only satisfy these conditions.
By the above presented algorithm we can prove the following theorem:
In other words there is at least one equilibrium point which is in
In the simulations of the above algorithm we can see the cell parameter projection on the (T, ∆, )-plane (Figure 6). We have red subregion in which we have three equilibrium points of our model and at least one is both stable and locally active; blue subregion in which we have either one or three equilibrium points and every equilibrium point is unstable; green subregion in which there is only one equilibrium point and it is both stable and locally active. By definition, red and green subregions in Figure 6a together constitute the edge of chaos. In Figure 6b we can see the plot of edge of chaos regime in the parameter (a, b, ) plane.
Through extensive numerical simulations we obtain that non uniform spatial patterns are generated in our CNN model with memristor synapses depending on initial conditions—see Figure 7.
Through the above numerical simulations, the following things were demonstrated: (a) excitable waves propagating on the memristor can modulate the memristor conductance which depends on the memristor’s polarity; (b) change of memristor conductance can modulate the velocity of the excitable wave propagation, and it is inversely proportional to the time constant of the model; (c) the model under consideration generates nonuniform spatial patterns which process depends on the initial condition of FitzHugh-Nagumo system (7) and (8), memristor polarity and stimulation.
3. Hysteresis CNN with memristor synapses
In this section we shall present mathematical study of hysteresis CNN (HCNN) with memristor synapses. In our model the cells are of first order and they have hysteresis switches. It is known from the literature [21, 22, 23, 24, 25, 26, 35, 37, 38, 39, 40, 41] that such models have many applications because they operate in two modes—bi-stable multi-vibrator mode and relaxation oscillator mode. We shall consider HCNN working in second one. When CNN operates in the relaxation oscillator mode then various patterns and nonlinear waves can be generated. Associative (static) and dynamic memories functions can be derived from the hysteresis CNN [35, 37, 41].
Let us consider hysteresis CNN with memristor synapses, which we shall call memristor hysteresis CNN (M-HCNN):
where denotes the state of the cell, the output is dynamic hysteresis function defined by:
, is defined as in which by we denote the memristance. When we insert memristor  in HCNN model we expect to obtain better resolution in static and dynamic images . We introduce a memristor in HCNN by replacing the original linear resistor. In this way it can exhibit nonlinear current-voltage characteristic with locally negative differential resistance. The main advantage of our memristor HCNN (M-HCNN) is the versatility and compactness due to the nonvolatile and programmable synapse circuits. In the circuit realization of M-HCNN the output function is not complex.
Let us consider M-HCNN model working in a relaxation oscillator mode described by
Below is the picture of relaxation oscillator under consideration (Figure 8).
3.1. Determination of edge of chaos domain in M-HCNN model
We shall apply theory of local activity [33, 34] in order to study the dynamics of M-HCNN model (15). The theory which will be presented below offers both constructive analytical and numerical method for obtaining local activity of M-HCNN. It is known [35, 41] that the cells of HCNN can exhibit complexity in the domain of cell parameters in which the cells are locally active. We shall develop constructive and explicit mathematical inequalities for identifying the region in the M-HCNN model (15) parameter space where complexity phenomena may emerge. By restricting the cell parameter space to the local activity domain we can achieve a major reduction in the computing time required by the parameter search algorithms. This will allow to determine and control chaos which will be useful for the future scientific and engineering applications [34, 41, 42].
(1) We chose the Laplace template of the following type in order to discretize the M-HCNN model (15). Then in relaxation mode the dynamics of an isolated cell when there are no control and threshold parameters can be written:
(2) We can find the equilibrium points of (16), which satisfy the equation . In general, this system may have four real roots as functions of the cell parameters.
(3) We calculate the cell coefficients , .
(4) We denote the trace of the Jacobian matrix at equilibrium point by and determinant by .
Remark. In order to provide physical implementation it is important to have appropriate circuit model for which we can use the results from the classical circuit theory. In order to obtain locally active cells it is sufficient the cell to act as a source of small signal in at least one equilibrium point. In this way the cell injects a net small signal average power into the passive resistive grids.
Let us now define stable and locally active region for the M-HCNN model (16).
Definition 1. We say that the cell is both stable and locally active region at the equilibrium pointfor M-HCNN model(16)if
This region in the parameter space is called.
(5) We shall define the EC region our M-HCNN model (16). According to [33, 34] it is such region in the cell parameter space where we can expect emergence of complex phenomena and information processing.
Definition 2. For M-HCNN model(16)edge of chaos region is such that there exists at least one equilibrium point both locally active and stable.
Based on the above algorithm we can prove the following theorem:
Theorem 2. We say that M-HCNN model(16)is working in edge of chaos regime if and only if the following conditions are satisfied: −1 < b < 3. In other words there is at least one equilibrium point which is locally active and stable.
We obtain the following edge of chaos domain for our M-HCNN model (16) through computer simulation:
The location of 16 cell parameter points arbitrarily chosen within the locally active domain. We have locally active and stable (or edge of chaos) in red; locally active and unstable in green; locally passive in blue (Figure 10).
3.2. Simulation results and some applications
In this section we shall consider two applications of M-HCNN model (16) in image processing. First one is for edge extraction and the second one is for noise removal. In our simulations we use programing environment MATLAB and we use a forward Euler algorithm with a time step size ∆t = 0.01. The dynamic hysteresis function h(x) can be programmed applying this algorithm is:
We shall start with the application of our M-HCNN model for edge extraction. The simulations are presented below (see Figure 11):
It is known that for feature extraction we firstly extract the edges of the image, which contain most of the information for the image shape. In the example provide on Figure 11 we show the original image—(a), and then the results which we obtain simulating M-HCNN model—(b) and standard CNN model—(c). It can be seen that the results from M-HCNN (16) model and CNN model are very similar.
Another application which we shall present is for noise removal. The results of our simulations are given on the figure below (see Figure 12):
The applications of CNN show that the linear image processing can be compared to spatial convolution with infinite impulse response kernels [21, 41, 42]. When taking the image by a camera from the real world there is a possibility it to be polluted with some noise. That is why noise removal is very important for CNN applications such as AI devices, IoT. In our case, the simulations of M-HCNN model (16) shown in Figure 12 present very good processing performance of noise removal similar to the simulations of standard CNN.
4. Conclusions and discussion
In this chapter, we stated the local activity theory for reaction-diffusion equations and hysteresis systems. However it can be generalized to other systems. In particular, the developed constructive procedure is applicable to any system whose cells and couplings are described by deterministic mathematical models. The crux of the problem is to derive testable necessary and sufficient conditions which guarantee that the system has a unique steady state solution at t → ∞. A homogeneous non conservative medium cannot exhibit complexity unless the cells, or the coupling network is locally active.
In the second part of the chapter we focus our attention on HCNN model which has memristor synapses. The concept of CNN is based on some aspects of neurobiology and is adapted to integrated circuits. CNN are defined as spatial arrangements of locally coupled dynamical systems, cells. The CNN dynamics is determined by a dynamic law of an isolated cell, by the coupling laws between the cell and by boundary and initial conditions. The dynamic law and the coupling laws of a cell are often combined and described by a nonlinear ordinary differential- or difference equation (ODE), the state equation of a cell. Thus a CNN is given by a system of coupled ODEs with a very compact representation in the case of translation invariant state equations. Despite of having a compact representation CNN can show very complex dynamics like chaotic behavior, self-organization, pattern formation or nonlinear oscillation and wave propagation. Analog CNN chip hardware implementations have been developed . The future of CNN implementation is in nano-structure computer architecture. CNN not only represent a new paradigm for complexity but also establish novel approaches to information processing by nonlinear complex systems. CNN have very impressive and promising applications in image processing and pattern recognition [22, 43]. After the introduction of the CNN paradigm, CNN Technology got a boost when the analogic cellular computer architecture, the CNN Universal Machine has been invented [24, 30]. The unique combined property of memristors  to store the information for a very long time after the power is switched off may allow the development and circuit implementation of memcomputing paradigms.
We develop algorithm for determination of EC domain of the cell parameter space for M-HCNN model (16). Two applications are presented—for edge detection and noise removal. The conclusions of the simulation results are that the image does not change when we vary the memristor weights which is possible because of the binary quantization of the output. The speed of the numerical simulations of our M-HCNN model could be enlarged due to the need of more iterations of the algorithm in order to obtain stable solutions. But the quality of the image does not change.
Acknowledgements are due to the project TE 257/25-1 between DFG and Bulgarian Academy of Sciences.
First author acknowledges as well the provided access to the e-infrastructure of the Centre for Advanced Computing and Data Processing, with the financial support by the Grant No BG05M2OP001-1.001-0003, financed by the Science and Education for Smart Growth Operational Program (2014–2020) and co-financed by the European Union through the European structural and Investment funds.