# System-Theoretic Methods for Designing Bio-Inspired Mem-Computing Memristor Cellular Nonlinear Networks

^{1}Fundamentals of Electrical Engineering, Institute of Circuits and Systems, Faculty of Electrical and Computer Engineering, Technische Universität, Dresden, Germany^{2}Department of Microelectronics, Brno University of Technology, Brno, Czech Republic^{3}Department of Electrical and Computer Engineering, University of California Santa Cruz, Santa Cruz, CA, United States^{4}Department of Electrical Engineering and Computer Sciences, University of California Berkeley, Berkeley, CA, United States

The introduction of nano-memristors in electronics may allow to boost the performance of integrated circuits beyond the Moore era, especially in view of their extraordinary capability to process and store data in the very same physical volume. However, recurring to nonlinear system theory is absolutely necessary for the development of a systematic approach to memristive circuit design. In fact, the application of linear system-theoretic techniques is not suitable to explore thoroughly the rich dynamics of resistance switching memories, and designing circuits without a comprehensive picture of the nonlinear behaviour of these devices may lead to the realization of technical systems failing to operate as desired. Converting traditional circuits to memristive equivalents may require the adaptation of classical methods from nonlinear system theory. This paper extends the theory of time- and space-invariant standard cellular nonlinear networks with first-order processing elements for the case where a single non-volatile memristor is inserted in parallel to the capacitor in each cell. A novel nonlinear system-theoretic method allows to draw a comprehensive picture of the dynamical phenomena emerging in the memristive mem-computing array, beautifully illustrated in the so-called Primary Mosaic for the class of uncoupled memristor cellular nonlinear networks. Employing this new analysis tool it is possible to elucidate, with the support of illustrative examples, how to design variability-tolerant bio-inspired cellular nonlinear networks with second-order memristive cells for the execution of computing tasks or of memory operations. The capability of the class of memristor cellular nonlinear networks under focus to store and process information locally, without the need to insert additional memory units in each cell, may allow to increase considerably the spatial resolution of state-of-the-art purely CMOS sensor-processor arrays. This is of great appeal for edge computing applications, especially since the Internet-of-Things industry is currently calling for the realization of miniaturized, lightweight, low-power, and high-speed mem-computers with sensing capability on board.

## 1 Introduction

On August 28^{th}, 2018 GlobalFoundries announced to halt the 7 nm chip development. After installing at least one Extreme-Ultraviolet Lithography (EUV) machine at one of its fabs, the foundry reckoned there would not be enough customers, interested in the cutting-edge 7 nm node technology process, to make chip development profitable (GlobalFoundries Ltd, 2018). Despite Intel, Samsung, and TSMC are still making efforts to reduce the size of *integrated circuits* (ICs) further, transistor scaling is approaching atomic boundaries, with an inevitable concurrent rise in manufacturing costs. This issue is known as *Moore wall*. With the Moore era (Moore, 1965) coming to a natural end (Williams, 2017), a great deal of resources have been deployed over the past decade toward the development of innovative disruptive nanotechnologies, which may enable the development of versatile multi-functional devices, that, opening the door toward the implementation of peculiar signal processing paradigms, would allow to boost the performance of conventional circuits and systems without the need to shrink the size of transistors any further. Two are the factors for the inefficiency of machines based upon the von Neumann architecture: 1) There is a large mismatch between processing time and shuttling time. This issue is known as *memory wall* or *von Neumann bottleneck*. 2) The energy dissipated by digital switching units is no longer following the exponentially decreasing trend, predicted by Landauer (Landauer, 1988), with the reduction in IC dimensions. This issue, known as *heat wall*, poses serious risks for the lifetime of transistors. Memristors (Chua, 1971; Chua and Kang, 1976) represent one of the most promising nanotechnologies to address the problems affecting state-of-the-art electronics. A current (voltage)-controlled non-volatile memristor (Chua, 2014; Chua, 2015) is a two-terminal device, whose resistance (conductance) can be tuned to some desired value by applying a current (voltage) signal through (across) it, and which remembers its resistance (conductance) after the current (voltage) source, in parallel to it, is disconnected (Chua, 2018a). It remembers its past! (Chua, 2018b). The most impressive and peculiar virtue of non-volatile memristors is the combined capability to store data, thanks to excellent data retention levels, achievable without the need for external batteries, and to process signals, through the rich nonlinear dynamics of the memory state, within a single nano-scale volume, which enables the implementation of *in-memory computing* and *mem-computing* paradigms^{1}, mimicking the distributed nature of memory and processing operations in the human brain, in future computing machines. Other distinctive qualities of memristors are low-power and high-speed operation, superior endurance, and, very importantly, good compatibility with CMOS technology. While in conventional memories data are stored as voltage levels, in memristors the physical quantity, which holds the information content, is the resistance. Given that all nonvolatile memories based upon resistance switching phenomena, irrespective of their constitutive materials and operating principles, are memristors (Chua, 2011), a wide range of nanotechnologies, including Resistive Random Access Memories (ReRAMs), Phase Change Memories (PCMs), Magnetic Tunnel Junctions (MTJs), Spin-Transfer-Torque Magneto-Resistive Random Access Memories (STTM-RRAMs), and Ferroelectric Tunnel Junctions (FTJs), are competing one with the other to produce the best performing data storage device for future brain-like computers. While many people believe that non-volatile nano-memristors will eventually replace conventional memories, including Flash Memories, Dynamic RAMs (DRAMs), and Hard Disk Drives (HDDs), the aforementioned nanotechnologies are not yet mature enough to draw conclusions on the portion of the nonvolatile memory market, which memristors will be able to cover in the next five years from now. However, *edge computing* technical systems already make use of memristive memories. Panasonic (Panasonic Ltd., 2013) have been launching mass production of micro-computers with 64 kB ReRAM-based data storage for battery-powered equipment, including portable devices for medical, healthcare, and security applications already in 2013, while Fujitsu (Fujitsu Ltd, 2019) has recently taken a step forward by offering 1024 kB ReRAMs for wearable units–e.g., smart watches and glasses–and hearing aids. Even when used simply as tunable resistors, memristors offer unique opportunities to enhance the performance of conventional data processing systems. Most computing tasks in artificial intelligence (AI) applications consist of machine-learning operations, such as object, image, and speech recognition, which require the calculation of a massive number of vector-matrix multiplications (VMMs). Nowadays these calculations are executed through expensive and bulky supercomputers. But with the advent of the memristor, which, leveraging Ohm’s law, naturally carries out a multiplication operation between the conductance it holds and the voltage falling between its terminals, outputting the result into the current flowing through it, it is possible to use a crossbar array (Li et al., 2018) to compute at unprecedented rates the product between a vector of voltages, distributed along the rows, and a matrix of memristor crosspoint conductances, with the computation result available in current form along the columns (refer to the Dot Product Engine (DPE) lab prototype developed at Hewlett Packard Enterprise (Hu et al., 2016)). Last but not least, given that the two constitutive elements of the human brain, namely the synapse and the neuron, are made of nonvolatile and volatile memristors, respectively, resistance switching memories allow to develop innovative neuromorphic circuits, which promise to outperform conventional purely CMOS counterparts in mimicking the functionalities of the human brain. Non-volatile memristor devices, in which the resistance may be finely tuned under excitation, may reproduce most closely the plastic response of biological synapses to external stimuli. Furthermore, there exist a large class of memristor physical nano-scale realisations, which, despite being unable to store data–for this reason they are classified as volatile memories –, feature the extraordinary capability to amplify infinitesimal fluctuations in energy (Bohaichuk et al., 2019; Kumar et al., 2020), a property which is referred to as *local activity* (Chua, 2005), and which enables the development of realistic electronic realisations of spiking neurons^{2} (Chua et al., 2012), the so-called neuristors.

Besides constitutive the ideal framework for modeling biological systems (Chua, 1998; Chua and Roska, 2002), *Cellular Nonlinear Networks* (CNNs) (Chua and Yang, 1988a; Chua and Yang, 1988b) represent a powerful multi-variate signal processing paradigm, which, featuring a bio-inspired architecture, operates in a massively parallel fashion, allowing to process data at very high rates, as necessary in time-critical *Internet-of-Things* (IoT) applications, nowadays. Purely CMOS analogue hardware implementations of the CNN signal processing paradigm are typically co-integrated with highly selective equal-sized sensor arrays to allow the solution of complex computing tasks directly where the acquisition of specific data takes place (Vázquez et al., 2018). A technological issue, which limits the applicability scope of these sensor-processor arrays, is related to the huge difference between the typically small minimum size of an element of the sensor matrix, and the relatively large minimum *integrated circuit* (IC) area, which a processing element of the CNN hardware realization usually occupies, due to the fact that it needs to accommodate memory units, which endow the resulting computing machine with local stored programmability on board, allowing to harness thoroughly the advantages associated with the massive parallelism of the CNN signal processing paradigm. The adoption of non-volatile memristors (Chua, 2015), capable to combine data processing and storage functionalities within a common nanoscale physical volume, in CNN circuit design may allow to increase significantly the spatial resolution of the cellular computing machine. Moreover, leveraging the rich nonlinear dynamics of resistance switching memories, the computing capabilities of the processing elements of a memristive CNN hardware implementation (Duan et al., 2015; Di Marco et al., 2017a; Di Marco et al., 2017b; Di Marco et al., 2018) may be extended beyond the operational boundaries of the cells of a traditional purely CMOS implementation.

CNNs process information through the analogue dynamics of the cells’ states, which converge toward distinct attractors depending upon inputs and/or initial conditions. While *wave-based computing*, where the cellular array carries out data processing tasks through the generation of specific dynamic patterns, is an active field of research (Weiher et al., 2019), there exists a huge library (Karacs et al., 2018) of image processing operations, which the nonlinear dynamic array may execute as the cells’ states approach predefined equilibria. This paper focuses on the performance of CNNs (M-CNNs) as *equilibria-based computing (mem-computing)* engine. Now, for a full exploration of the potential of memristors in electronics, recurring to concepts from *nonlinear system theory* is necessary. In fact, linear system-theoretic methods are not suitable for the analysis and design of memristor-based circuits. However, as is the case here, converting traditional nonlinear circuits to memristive equivalents may require the extension of classical nonlinear system-theoretic techniques. The *Memristor Cellular Nonlinear Network* (M-CNN), proposed in Tetzlaff et al. (2020), differs from a standard time- and space-invariant two-dimensional CNN (Chua and Yang, 1988a; Chua and Yang, 1988b), characterized by first-order cells, and typically implemented in hardware (Vázquez et al., 2018), for the inclusion of a single non-volatile memristor in parallel to the capacitor in the circuit implementation of each processing element. One of the most powerful tools for the analysis of nonlinear dynamical systems with one degree of freedom is the *Dynamic Route Map* (DRM) (Chua, 2018a), which represents the system-theoretic technique of reference for the investigation of CNNs with first-order processing elements. Since the memristive cell in the proposed M-CNN features two degrees of freedom, the investigation of the cellular array calls for the generalization of the DRM graphical tool, applicable to first-order systems only. The modified DRM graphical tool, applicable to second-order dynamical systems, is known as *Second-Order Dynamic Route Map* (DRM_{2}) (Tetzlaff et al., 2020). The application of this novel system-theoretic technique to the model of the proposed M-CNN allows to gain a deep insight into the rich nonlinear behaviour of its second-order processing elements, unveiling dynamical phenomena, which may not emerge in the original cellular array (Ascoli et al., 2020b). The DRM_{2} graphical tool lies at the basis of a systematic methodology to design variability-tolerant mem-computing arrays with second-order memristive cells (Ascoli et al., 2020a).

The structure of the paper is organized as follows. Section 2 revisits the theory of CNNs, explaining the invaluable role of the classical DRM graphical technique to analyze standard arrays of locally coupled processing elements, and elucidating through an illustrative example the traditional method, based upon this system-theoretic tool, to program the cellular computing engine for the execution of a predefined image processing task. Section 3 first defines the class of M-CNNs under study, including the model of the non-volatile memristor hosted in each cell, secondly extends the DRM graphical tool to second-order dynamical systems, elucidates how this allows to draw a comprehensive picture of the nonlinear dynamics of each memristive processing element, and finally presents a rigorous procedure, based upon the DRM_{2} system-theoretic technique, to design cellular mem-computing structures with second-order memristive cells. Sections 4 and 5 are devoted to the application of the M-CNN design methodology for operating the multifunctional memristive cellular computing engine as image processing system and as memory bank, respectively. A brief discussion, summarizing the significance of the research work, is provided in section 6. Conclusions are finally drafted in section 7.

## 2 Analysis and Design of Memristor Cellular Nonlinear Networks

Cellular Nonlinear Networks (CNNs) constitute a bio-inspired multivariate signal processing paradigm, which, based upon a massively parallel information flow, enables computations at very high rates, and is amenable to a Very Large Scale Integration (VLSI) circuit realization, which, centered around a non-von-Neumann machine architecture, enables computational universality. Introducing memristive devices in CNN VLSI design may provide two main benefits. Firstly, the rich spectrum of nonlinear dynamic phenomena, appearing in resistance switching memories, may simplify or extend the functionalities of traditional CNNs. Secondly, the unique combined capability of nonvolatile memristors to compute and store data within the same nanoscale physical medium may render unnecessary the need to include spacious data storage units within the circuit implementation of each cell, allowing to improve considerably the number of processing elements fitting into the IC design area allocated to the non-von-Neumann computing machine.

### 2.1 Theory of Cellular Nonlinear Networks

The theoretical foundations of CNNs were laid in 1988 by L. Chua (Chua and Yang, 1988a; Chua and Yang, 1988b). In the most general case, a CNN consists of a spatially discrete collection of locally coupled *k*^{th}-order continuous-time processing elements, called cells, arranged at regular positions within a *l*-dimensional lattice. The architecture of a small two-dimensional CNN with *M* = 6 rows and *N* = 6 columns is presented in Figure 1A, under the assumption that each cell ^{3}. Each cell is assigned a state, an input, an output, as well as a threshold. The rate of change of the state of a cell is influenced by the inputs and outputs of its 8 adjacent neighbors, as well as by its own input and output, as respectively sketched through eight brown directed segments and through one magenta directed loop in Figure 1A for the processing element located where the 3^{rd} row crosses the 4^{th} column. The block diagram in plot (b) of Figure 1 illustrates once more the key factors affecting the dynamical behaviour of a CNN cell. The neighbors’ inputs (outputs) are modulated by feedforward (feedback) synaptic weights before accessing the cell to affect the time evolution of its state, similarly as it occurs in biological neural networks. The cell ^{4} (

in case it is physically coupled to its 8 adjacent neighbors only^{5} i.e., *r* = 1. With reference to the circuit of Figure 2, the main physical quantity within the input stage, the computing core, and the output stage respectively are the input voltage

where

generally known as *standard nonlinearity*, where *negative saturation region**linear region**positive saturation region*

where *z* is a dimensionless parameter referred to as *threshold* in CNN theory (Chua and Roska, 2002), while *I* is a coefficient with positive 1 A value, symbolizes the threshold current. Further, ^{6}, while, importantly,

where *feedback (feedforward) synaptic weight* in CNN theory (Chua and Roska, 2002), represents the feedback (feedforward) synaptic current due to the neighboring cell

**FIGURE 1**. **(A)** Physical connectivity among the locally coupled cells of a two-dimensional CNN with six rows and six columns. **(B)** Schematic illustration of the main features of a CNN cell, revealing some of its analogies with a biological neuron, which explains the nomenclature originally introduced to characterize the locally coupled nonlinear dynamic array, namely *Cellular Neural Network*.

**FIGURE 2**. Input stage, computing core, implementing the state Eq. 1, and output stage of the circuit realization of a standard space-invariant CNN cell

Typically, the right hand side of the CNN cell state Eq. 1 may be recast as

where *Internal Driving Point (DP) Component* is a function of the cell state, being defined as

while *offset current*, mostly accounts for the coupling effects, being expressed by

where *z*, the feedback synaptic weights *gene* in CNN theory. A gene defines the set of rules, which apply concurrently in the neighborhood of each cell, allowing the standard space-invariant CNN to carry out a given computation.

Remark 1. A CNN may be used to carry out any computing task. The calculations are based upon the analogue dynamics of the cell states. As transients vanish, depending on cell parameter settings, inputs and initial conditions, each capacitor voltage may tend toward an equilibrium, or converge to an oscillatory waveform, which may be periodic, quasi-periodic, or even chaotic. A CNN may then process the information, inserted as cell inputs and/or initial conditions, in two different forms i.e., producing computation results in the form of equilibria or waves, hence the names CNN equilibria-based computing or CNN wave computing attributed to the respective operating principle. In this manuscript the attention is focused on CNNs computing via equilibria^{7}*.* Let us elucidate the classical method to design a CNN, so that it may execute a fundamental image processing task^{8}, carrying out the result of the computation in the cell equilibria. A rigorous technique to synthesize the gene of a standard CNN, so as to allow the successful execution of a given equilibria-based computing task, even in the presence of deviations of some cell circuit parameters from their nominal values, was introduced in (Zarándy, 2003), and comprehensively presented in (Itoh and Chua, 2003). Before summarizing the line-of-thought at the basis of this classical methodology, let us provide a brief overview of the nonlinear dynamics of a CNN processing element.

### 2.2 Key Features of the CNN Cell Dynamics

The first aspect to consider for the synthesis of a suitable CNN gene is the choice of a proper value for the self-feedback synaptic weight *Internal DP Characteristic*, consisting of the arrowed locus of the rate of change of the state ^{−1}, arrows, supeimposed over the Internal DP Characteristic, point toward the east (west) to indicate an increase (a decrease) in the state ^{9} Figure 3, plot (a) shows how

• If

• If

• If

**FIGURE 3**. **(A)** Family of Internal DP Characteristics, which a CNN cell admits for each value of the self-feedback synaptic weight **(B)** Family of Shifted DP Characteristics, which a CNN cell admits for

The ordinates of the two breakpoints of the three-segment^{10} piecewise-linear Internal DP Characteristic, located at *Shifted*^{11}*DP Characteristic* (Chua and Roska, 2002), i.e. the locus of the state rate of change ^{12} CNN is *completely stable* (Chua and Roska, 2002), in the sense that the state *bistability criterion* (Chua and Roska, 2002), in case the condition

holds true, the slope of the Internal DP Characteristic in the linear region–refer to Eq. 9 – is positive, and, consequently, the stable equilibria of each cell under *steady state* as far as its cell output voltages are concerned^{13}. It follows that, irrespective of the location of the CNN cell

With reference to the viewgraphs of Figure 3B, neglecting the marginal cases, three are the possible equilibria configurations, which a cell

• If

• If

• If

### 2.3 A Systematic DRM-Based Methodology to Design Robust CNNs

The standard method (Zarándy, 2003; Itoh and Chua, 2003) to synthesize a suitable CNN gene for the execution of a given computing task is based upon the set up and later solution of an ad-hoc set of inequalities, expressed in terms of unknown gene parameters, ensuring a robust accomplishment of the computing task of interest. The functionalities of a given uncoupled^{14} standard space-invariant completely stable CNN, satisfying the hypotheses of the bistability criterion, are dictated by a set of local rules, establishing the asymptotic value for the state

For the reasons motivated above, in CNN designs for image processing applications, it is common to set the numerical value for the self-feedback synaptic weight ^{15}, and, consequently, the emergence of CNN computation errors. The next step is to write down an inequality, establishing a constraint for the offset current, and ensuring that, in such critical scenario,

### 2.4 Edge CNN

The aim of this section is to synthesize the gene of a standard space-invariant non-autonomous uncoupled CNN so that the resulting nonlinear dynamic array is able to extract the edges from an input binary image. Next, the classical CNN design methodology, briefly described earlier, will be applied to the cell ODE (1) in order to achieve this purpose. The three local rules, which each cell should obey, are reported in Table 1, where *i* and column *j* in the input binary image under rule 1. Here 3 of the 8 neighboring pixels are black. The pixel in the position *i* and column *j* in the input binary image under rule 2 (3). Here all (4) of the 8 neighboring pixels are black. The pixel in the position

**TABLE 1**. Local rule triplet, which, irrespective of its location within the cellular array, a processing element

**FIGURE 4**. Graphical illustration of the application of the EDGE CNN local rules 1 for **(A)**, 2 **(B)**, and 3 for **(C)**. Each of the three plots visualizes, on top, a *i* and column *j* in the input binary image, and, below, the pixel at position

As clarified by Figure 5A, the CNN under design is expected to extract the edges from an input binary image, visualizing them in the output binary image at steady state. Since the CNN is meant to be uncoupled, the offset current from Eq. 11 reduces to^{16}

**FIGURE 5**. **(A)** Graphical illustration of the operating principles of the CNN under design. **(B)** EDGE CNN SDR for zero offset current. Here *b* value for each of the feedforward synaptic weights, except for

Assuming that all the feedforward synaptic weights, with the exclusion of

where the argument reveals the dependence of the offset current upon *z*. Figure 5B shows the directed Internal DP Characteristic^{17}. The state Eq. 1 under

Let us set the initial condition

In order to fulfill rule 1, where

It is worth pointing out that the farther away from the horizontal axis, within the plane lower half, would be positioned the right breakpoint of the ^{18}.

**FIGURE 6**. Graphs clarifying the line of reasoning behind the DRM synthesis strategy adopted in Itoh and Chua (2003) to select a suitable gene allowing the resulting CNN to apply the local rule triplet of the binary image edge extraction operation in the 9-cell neighborhood of each processing element. The worst-case scenario in rule 1 is analyzed in **(A)**, where **(A)** allows to investigate rule 2 as well. The worst-case scenario in rule 3 is illustrated in plot **(B)**, where **(A)**, in the worst-case scenario from rule 1 (in rule 2) the cell state **(B)**, in the worst-case scenario from rule 3,

Let us now derive the condition allowing the CNN to apply rule 2 from Table 1 in the sphere of influence of any processing cell

In order for the CNN under design to apply rule 3 from Table 1 in the ^{19}

For a robust CNN design the right breakpoint of the ^{20}.

For the parameter setting reported in the caption of Figure 5, the three inequalities Eqs. 15–17 are solved through a geometric approach on the *z*–*z* and

**FIGURE 7**. Illustration of the geometrical analysis adopted to solve the three inequalities Eqs. 15–17, derived through the classical CNN design method (Itoh and Chua, 2003) to synthesize a suitable gene for a cellular array, intended to extract edges from a given input binary image, and reducing to

### 2.5 Limitations of the CNN Paradigm and of Its Hardware Implementation

Since each of their processing elements interacts simultaneously with the respective neighbors, CNNs may process multi-variate signals in a massively parallel fashion, as crucially necessary in time-critical application fields, such as industry process control, electronic surveillance, medical augmented reality, and IoT smart sensing. In order to harness more efficiently the bio-inspired operating principles of these nonlinear dynamic arrays, which make them a suitable mathematical framework for modeling neural systems, Chua and Roska proposed an innovative computer, called CNN Universal Machine (CNN-UM) (Roska, 1993), to implement their signal processing paradigm. The CNN-UM, fabricated in various forms over the years through the well-established CMOS technology^{21} (Vázquez et al., 2018), consists of an array of locally coupled computing units, each of which is endowed with data storage blocks, which allow to distribute the memory across the cellular array, endowing the computing machine with a truly non-von Neumann architecture, and to reconfigure the array so as to solve any computation problem. Thanks to their massively parallel computing power, CNN-UM hardware realizations (Vázquez et al., 2018) may process images at rates as high as 30,000 frames per second. Considering that, furthermore, a universal cellular array may be physically realized within the IC area of a single chip (Vázquez et al., 2018), CNNs are particularly suitable for the development of miniaturized IoT technical systems, in which the integration between a matrix of sensing elements and a network of locally coupled computing units with local stored programmability on board enables information processing at the same location, where data detection takes place. A major problem, which prevents to widen the applicability scope of this class of sensor-processor arrays, is the limited degree of complexity of the dynamical phenomena, which may possibly emerge within their physical media, due to the simplicity of the input-output behaviours of the electrical components employed in the CNN-UM constitutive blocks. Thanks to their extremely rich dynamics, memristors may be adopted in novel designs of cellular computing arrays so as to extend significantly the spectrum of asymptotic spatio-temporal behaviours, which purely CMOS CNNs may currently exhibit. Another critical issue, which affects the performance of technical systems, combining sensing and processing functionalities on the same physical platform, is due to the rather low spatial resolution of state-of-the-art CNN-UM hardware realizations, originating from the presence of spacious data storage units within their computing units, as discussed earlier. This limits the maximum number of sensing and processing elements, which may be paired^{22} one-to-one within the available IC area of these IoT commercial products (Toshiba Ltd., 2012). The adoption of memristive devices, endowed with memprocessing capabilities, may allow to obviate the inclusion of additional memory banks within the IC design of each CNN-UM computing unit, allowing to shrink considerably its size, and enabling the future realization of sensor-processor arrays with unprecedented spatial resolution, of great appeal to the IoT industry, nowadays. In this respect, it is timely to commence investigations aimed to explore the functionalities of Memristor CNNs (M-CNNs). In general, introducing memristors in the circuit implementation of a CNN processing element^{23} increases the order of its ODE model, calling for the development of a new theory to investigate the operating principles of the resulting nonlinear dynamic array, and to program its gene to allow the accomplishment of a pre-defined memcomputing task. The theoretical foundations of M-CNNs shall be discussed in the section to follow.

## 3 Theory of Memristor Cellular Nonlinear Networks

Memristors are the key technology enabler for the hardware implementation of innovative memcomputing paradigms. This section provides some evidence for this claim, establishing the theoretical foundations of a class of cellular memprocessing structures, which we call M-CNNs, as anticipated in section 2.5. In order to realize one of the proposed M-CNNs a first-order non-volatile memristor^{24}

**FIGURE 8**. Circuit implementation of the M-CNN cell

The next section reports the mathematical description of the proposed memristive cellular array.

### 3.1 M-CNN Model

The M-CNN cell ^{25} (

The first ODE Eq. 18 governs the time evolution of the state *generic memristors*, defined by the DAE set

Note that, within the processing element

The state evolution function

The memristor state

where

**FIGURE 9**. **(A)** Course of ^{−1} s^{−1}, ^{−1} s^{−1}, and **(B**,**C)** Window functions, appearing in the state evolution function Eq. 23, and preventing **(upper)** bound

Since the memristor state existence domain ^{26} in the Pershin and Di Ventra model (Pershin et al., 2009), adopting continuous and differentiable functions, inspired to Biolek’s window (Biolek et al., 2009), and reading as

where

in which Eqs. 22, 24 were employed to model the cell memristor current *phase portrait* concept from the theory of nonlinear dynamics (Strogatz, 2000), the next section introduces a new system-theoretic notion, which we name *Second-Order DRM* (DRM_{2}), enabling the investigation of the memcomputing capabilities of cellular nonlinear arrays with second-order memristive cells.

**FIGURE 10**. **(A)** SDRs foliating from the memristor DRM for **(B)** Proof of evidence for the analogue dynamic response of the memristor state, hosted by the cell **(C)** Pinched hysteresis loop emerging on the **(B)**. The memristor model parameters are set as follows: ^{−1} s^{−1}, ^{−1} s^{−1},

### 3.2 A Generalized DRM Technique for the Analysis of M-CNNs With Second-Order Processing Elements

In this section we extend the classical DRM methodology (Chua, 2018a) for the analysis of a nonlinear dynamic system with two degrees of freedom. Focusing, in particular, on the second-order M-CNN cell under study, the *State Dynamic Portrait* (SDP).

Remark 2. With reference to the qualitative drawing in Figure 11, a SDP is a two-dimensional graph associated to a prescribed choice for the offset current value. It may be obtained as follows. First, the phase plane *–*

• Green region I:

• Yellow region II:

• Cyan region III:

• Gray region IV:

Then the loci

• Red crosses:

• Magenta diamonds:

• Black circles:

Particularly, the local instability (stability) of an equilibrium, studied by linearizing the state equations and studying the properties of the Jacobian, is graphically illustrated in a given SDP by means of a hollow (filled) black circle. The dynamical behaviour of the state variables from any initial condition of interest may be qualitatively inferred by inspecting the direction of the vector field ^{27}, moving through regions I, II, III, and IV, proceed in the south-west, north-east, north-west, and south-east directions, as time goes by, respectively. The numerical integration of the pair of first-order coupled ODEs Eqs. 18, 19, for initial conditions in the set of interest, allows to confirm this qualitative investigation on a quantitative basis, allowing to endow the partitioned plane, already accomodating nullclines and equilibria, with a number of phase plane trajectories, extracted by plotting the two solutions _{2}). Finally, the family of SDPs, obtained for each offset current value within a certain set of interest, takes the name of Second-Order DRM (DRM_{2}). The proposed generalized DRM methodology may be used to analyze the operating principles of a given M-CNN with second-order memristive cells. Most importantly, the DRM_{2} graphical tool allows to develop a systematic procedure to program one of the memristive cellular arrays under focus for the execution of a predefined memcomputing task, as outlined in the next section.

**FIGURE 11**. Exemplifying SDP, which the M-CNN cell

Remark 3. The DRM_{2} graphic tool features a much more general applicability scope than this paper demonstrates. In fact, it allows to investigate any second-order dynamical system, including memristive circuit elements with two degrees of freedom.

### 3.3 A Rigorous DRM_{2}-Based Methodology for Robust M-CNN Design

The proposed DRM_{2}-based M-CNN design methodology (Ascoli et al., 2020a) allows to program the memristive nonlinear dynamic array i.e., to choose numerical values for the 19 cell core parameters^{28}^{29} (Chua, 1998), which, depending upon the specific data storage or processing operation to be executed, dictate the steady-state value^{30} of its output voltage

1. On the basis of the memcomputing task assigned to a given M-CNN, and with reference to the processing element _{2}, under target.

2. In order to derive numerical values for the parameter set _{2}, a number of inequalities, constraining, for each scenario of any rule, the behaviour of the ^{31} through the use of the second-order ODE system Eqs. 18, 19, with the expression for the offset current

3. A set of cell parameter values, satisfying concurrently all the aforementioned inequalities, shall be determined by means of a graphical approach, or through a numerical algorithm, depending upon the number of unknowns. Integrating numerically the state Eqs. 18, 19 of the M-CNN cell

### 3.4 Application of the M-CNN Design Methodology to Execute Fundamental Memcomputing Tasks

In this section the proposed cell DRM_{2} synthesis technique is applied to the cell model Eqs. 18, 19 to program the M-CNN to execute an image processing operation and a couple of memory functions, namely the data storage and retrieval tasks. Before presenting the M-CNN design examples, it is instructive to identify the most important properties of the second-order system Eqs. 18, 19 through the application of fundamental concepts from the theory of nonlinear dynamics (Strogatz, 2000).

From the first M-CNN cell ODE Eq. 18, the formulas for the

Employing now the second M-CNN cell ODE Eq. 19, the

for

for

for

Remark 4. As it follows from Eq. 29, under

also of the two disjoint sets

**FIGURE 12**. Cell SDP, emerging for the fixed circuit parameter setting from Table 1, under **(A)**, a stable isolated equilibrium, as well as a line of equilibria with stable (unstable) character to the left **(right)** of a bifurcation point for **(B)**, and, finally, two stable isolated equilibria, as well as a continuum of unstable equilibria for **(C)**.

Remark 5. The application of the proposed design method to the specific M-CNN cell model Eqs. 18, 19 is unable to control existence and/or massage the shape of the _{2} for the accomplishment of a predefined memcomputing task. The equilibria, lying at the intersections between the

if

in the negative saturation region, at

if

as well as at

if

in the linear region, and at

if

in the positive saturation region.

Remark 6. Under

represents a possible equilibrium for the M-CNN cell in the linear region. Moreover, in case

and decreases if

Thus, as revealed by the illustrative cell SDP example of Figure 11, the motion of a trajectory point

for

for

for

##### 3.4.1 Zero Offset Current Scenario

It may be proved that, unlike a standard CNN processing element, the M-CNN cell *monostability* for

**TABLE 2**. Location and local stability property of each of the equilibrium points, which a M-CNN cell may possibly admit, depending upon

Specifying the values^{32}, reported in Table 3, for all the fixed parameters of the cell circuit of Figure 8, the viewgraphs in plots (a), (b), and (c) of Figure 12 illustrate the SDP of a M-CNN processing element, accommodating no linear resistor in the memcomputing core, under zero offset current, and for a specific value of the self-feedback synaptic weight

**TABLE 3**. Setting of specific M-CNN cell circuit parameters, specifically α, β, and *p* from Eqs. 26, 27, *I* from Eq. 4,

##### 3.4.2 Non-Zero Offset Current Scenario

Allowing a non-null offset current, accounting mostly for the coupling effects, to flow through the capacitor in the circuit of Figure 8 may endow the processing element with monostability, which is useful for the accomplishment of certain mem-computing tasks, as will be clear from the discussion of some M-CNN designs in the sections to follow. Table 4, in which

classifies the number, location, and local stability property of all the equilibria which a M-CNN cell may possibly admit for all the possible combinations of self-feedback synaptic weight

**TABLE 4**. Location and local stability property of each of the equilibrium points, which a M-CNN cell may possibly admit, depending upon

Remark 7. Interestingly, this table allows to draw the codimension-2 bifurcation diagram of Figure 13, in which, without loss of generality, *was set to**.* This graph, which, taking inspiration from CNN theory (Chua, 1998; Chua and Roska, 2002) is called M-CNN Primary Mosaic, visualizes the partitioning of the *–**.*In this manuscript the proposed DRM_{2}-based M-CNN design methodology shall be applied to the model Eqs. 18, 19 of a cell belonging to the class of uncoupled M-CNNs, and featuring an offset current, which, in comparison to its most general formula, namely Eq. 11, reduces to Eq. 13.

**FIGURE 13**. The *M-CNN Primary Mosaic*: codimension-2 bifurcation diagram illustrating all the admissible equilibria the two states of the second-order cell from the class of uncoupled M-CNNs may approach asymptotically depending upon the specific region of the

## 4 M-CNN as a Bio-Inspired Image Processing Engine

A M-CNN may be programmed to carry out any image processing operation, which a classical CNN is able to execute. To provide some evidence for this claim, the next section discusses the system-theoretic design of a memristive cellular array for the extraction of edges from a binary image.

### 4.1 Edge M-CNN

This section is devoted to the design of a ^{33}, are reported in Table 1 from section 2.4. In order to ensure that the memprocessing elements obey this local rule set, it is wise to synthesize the cell SDP pertaining to each scenario from rules 1 and 2 (rule 3) in such a way that it accommodates one and only one equilibrium located in the negative (positive) saturation region i.e., ^{34} rule 1 for _{2} synthesis strategy, in any scenario of rule 1 and for rule 2 (under all circumstances in rule 3), a state vector _{2}, as illustrated in Figure 15. Previous to initiate this investigation, a couple of aspects should be pinpointed. Firstly, the cell ODE initial condition

**FIGURE 14**. **(A)** SDP of a cell **(B**,**C)** SDP of a cell

**FIGURE 15**. Qualitative visualization of the strategy adopted to massage the shape of the cell DRM_{2} in the EDGE as well as in the STORE M-CNN designs. Here the stepwise application of the proposed systematic gene synthesis methodology in the linear region (a.1) ((b.1)), in the negative (positive) saturation region (a.2) ((b.2)), and in the positive (negative) saturation region (a.3) ((b.3)) of the phase plane ^{−1} loci, depending upon the sign of

##### 4.1.1 Edge M-CNN Cell DRM_{2} Synthesis in the Linear Region

With reference to plot (a.1) (b.1) of Figure 15, the aim of this section is to make sure that, under all circumstances in rule 1 and for rule 2 (in all scenarios of rule 3), the characteristic

ensures a positive sign for the denominator of the rational function on the right hand side of Eq. 35 irrespective of the value assumed by the memristor state ^{35} in each scenario of rule 1 and for rule 2 (under all circumstances in rule 3) via

ensures that the

Let us now study the direction of motion of the state vector

Since, as established by constraint Eq. 57, in rules 1 and 2 (rule 3) the right hand side of inequality Eq. 58 assumes values larger (lower) than ^{36} the gray (yellow) region IV (II) or the green (cyan) region I (III), respectively, as indicated in plot (a.1) ((b.1)) of Figure 15.

##### 4.1.2 Edge M-CNN Cell DRM_{2} Synthesis in the Saturation Regions

The goal of this section is twofold. On one hand, in each scenario of rule 1 and for rule 2 (under all circumstances from rule 3) the cell SDP is expected to accommodate one and only one GAS equilibrium point, specifically Eqs. 39, 45, over the domain

in each of the first (latter) set of scenarios, ensures that this unique ^{−1} locus lies within the negative (positive) saturation region for all *b* of the off-center synaptic weigths in the feedforward template, it may be easily realized that, under *z*, and *b*, respectively set to ^{37} in parallel to the capacitor in the memcomputing core of Figure 8, the two conditions, respectively descending from constraint Eqs. 56, 57 under the first inequality sign option for the worst-case scenario _{2} in the worst-case scenario ^{38}, upon setting the resistance of each memristor to

**FIGURE 16**. **(A)** Geometry-based approach to the solution of the system of three non-redundant inequalities, resulting from the application of the systematic M-CNN design methodology from section 3.3 for the accomplishment of the binary image edge extraction task. The magenta region contains the set of admissible solutions of the constraint triplet, which consist of inequality Eq. 55, and of the couple of conditions descending from the constraint pair Eq. 57 under the first inequality sign option in the worst-case scenario *z*, *b*, and **(A–C)** of Figure 13 show the SDPs foliating from the resulting EDGE M-CNN cell DRM_{2} in the worst-case scenario **(B**–**D)** Proof of evidence for the proper functionality of a **(D)** depicts the steady-state output binary image of the EDGE M-CNN, once the black-and-white image, illustrated in plot **(B)**, is loaded to its input, a zero is assigned to the voltage falling across each capacitor at the beginning of the simulation, as visualized through the gray image in plot **(C)**, the initial condition on the resistance of each memristor is set to 5 k

Remark 8. The insertion of a single memristor within the circuit implementation of the cell of a standard time- and space-invariant CNN allows to endow the resulting memristive array with novel functionalities, including the capability to read and write data locally within each processing element without the need to accommodate additional memory units, which are currently responsible for the poor spatial resolution of state-of-the-art CNN-UM hardware realizations.

## 5 M-CNN as a Memory Bank: Write/Read Functionalities

The operating principles of a M-CNN programmed to write or read input binary data into or from the resistances of its memristors are elucidated below.

### 5.1 Store M-CNN

The aim of this section is to synthesize the gene for programming the cell ^{39}, for all

**TABLE 5**. Pair of local rules, which the M-CNN cell

Taking inspiration from the strategy adopted earlier on in the synthesis of a suitable gene for programming the bio-inspired memristive array to extract edges from an input binary image, a possible approach to design the STORE M-CNN is to make sure that the cell SDP accommodates one and only one globally asymptotically stable equilibrium i.e.,

In order to facilitate the comprehension of the STORE M-CNN design, let us anticipate its outcome. Plot (a) ((b)) in Figure 17 depicts the cell SDP synthesized through the proposed DRM_{2}-based system-theoretic technique to store the input voltage ^{−1} locus Eqs. 34, 36, indicated via magenta diamonds and located in the negative (positive) saturation region, and that the ^{−1} locus Eqs. 31, 33, identified via red crosses in the phase-plane lower (upper) half, in the equilibrium Eqs. 39, 45, where, as marked by means of a black circle, the memristor stores the off (on) resistive state _{2} this way, any trajectory of either SDP, visiting the region below (above) the single ^{40} non-autonomous array is expected to suffice for the accomplishment of the binary data writing task, the expression for the offset current in Eq. 13 reduces to

Under the hypothesis that a value is preliminarily assigned to the self-feedback synaptic weight _{2} system-theoretic method will aim to the determination and later solution of an ad-hoc IS in the pair of unknown parameters *z*. The proposed gene synthesis technique is first applied to derive the necessary constraints for well-behaved phase-plane trajectories in the linear region of each cell SDP.

**FIGURE 17**. **(A**,**B)** SDP of the processing element *i* and column *j* of a given input binary image, in the form of off (on) resistance *D*, as expected from rule 1 (2) of Table 5. **(C)** Graphical determination of the solutions of the non-redundant inequality pair Eq. 62, descending from the application of the proposed DRM_{2}-centered system-theoretic bio-inspired array design method to the model Eqs. 18, 19 of the processing element *z*–_{2}, foliating in the SDPs of plots **(A)** and **(B)**, associated to a nonzero offset current **(C)**, and located well within the magenta region, conveniently away from its boundaries. For the sake of completeness, with reference to Figure 14 (a.3) ((b.3)), the final solution of the non-redundant inequality system defines the shape of the characteristic

##### 5.1.1 Store M-CNN Cell DRM_{2} Synthesis in the Linear Region

The purpose of the following mathematical derivations is to make sure that, under negative (positive) one V-valued cell input voltage, the locus of points, lying on the characteristic

the graph of the function

Under rule 1 (2), the whole phase-plane region _{2} in the phase-plane saturation regions.

##### 5.1.2 Store M-CNN Cell DRM_{2} Synthesis in the Saturation Regions

Two are the aims of the mathematical treatment to follow. Firstly, referring to Figure 15 (a.2) ((b.2)), we shall make sure that, under ^{−1} locus, which the cell SDP may ever accommodate under rule 1 (2), we shall impose that the characteristic ^{−1} locus Eqs. 34, 36 exhibits upward (downward) concavity in its monotonic decrease (increase) with the memristor state. Thus, making sure it lies below (above) the horizontal line

under

Basic mathematical analysis reveals that, under the hypothesis of rule 1 (2), the function ^{41}*z*–*z*, and

Setting the values for the core circuit parameters of each cell of a M-CNN, featuring *i* and column *j* of a binary image, featuring the same spatial resolution as the M-CNN, and shown in Figure 18A, is first mapped onto a negative (positive) one V-valued input voltage

**FIGURE 18**. **(A)** Input binary image with **(B**,**C)** Black-and-white picture coding the random initial condition assigned to the capacitor voltage (memristor state) in each cell of a M-CNN with same rows and columns as the input image. **(D)** Graphical illustration of the data stored in the memristive memory at the end of the data writing operation. The binary picture in **(D)** appears to be the logically inverted version of the black-and-white image in **(A)**, due to the convention, intentionally adopted here, to map each white (black) pixel of the input binary image to a negative (positive) one V-valued input voltage for the corresponding M-CNN cell, and to code the memory content of each off (on) memristor at equilibrium through a black (white) pixel.

### 5.2 Recall M-CNN

The purpose of this section is to shape the DRM_{2} of the processing element ^{42} output voltage

**TABLE 6**. Set of local rules, which are imposed on the processing element

Given that the RECALL M-CNN is autonomous, differently from the strategy adopted in the STORE M-CNN design, the gene synthesis approach, followed in this section, aims to massage one and only one SDP, hosting the solutions of the second-order ODE Eqs. 18, 19 for both local rules from Table 6, and constituting, as a result, the DRM_{2} itself. In order to achieve this purpose, the cell should be programmed so as to operate as a bistable dynamical system: in case the memristor, it accommodates, sits in the off (on) resistive state, the vector field flow should guide the trajectory toward the equilibrium

**FIGURE 19**. **(A)** SDP of the bistable cell **(B)** Graphical illustration of the geometric analysis, carried out in the parameter plane *z*–_{2} of each M-CNN processing element in such a way to retrieve the information stored in the memristor **(A)**, was derived for the particular solution **(B)**, resides, to some extent, away from the white parameter plane region, where at least one of the three conditions does not hold true.

With reference to the cell SDP in Figure 19A, in our strategy we first imposed the existence of a _{2} in the linear, negative saturation, and positive saturation region, as desired (refer, once again to Figure 19A), are visualized through illustrative viewgraphs in Figures 20A–C, respectively.

**FIGURE 20**. Qualitative sketches illustrating graphically the steps of the system-theoretic method to massage the cell SDP in the linear **(A)**, negative saturation **(B)**, and positive saturation **(C)** regions for the derivation of a suitable gene to program the M-CNN for the data retrieval task.

A rigorous mathematical analysis of the cell ODE Eqs. 18, 19 allows the derivation of a suitable IS for the creation of an ad-hoc cell SDP, combining the coloured phase-plane regions in plots (a), (b), and (c) of Figure 20. A preliminary requirement for initiating the investigations is to fix the expression for the offset current. Conjecturing that the use of a simple isolated autonomous M-CNN would be sufficient for the accomplishment of the data recall operation, simplifying Eq. 13, the following formula may be assigned to

Assuming that the conductance *z*, which will enable the determination of the domain of admissible solutions on the basis of a geometrical analysis. Let us commence the systematic mathematical treatment from the linear region of the standard nonlinearity of Eq. 3.

##### 5.2.1 Recall M-CNN Cell DRM_{2} Synthesis in the Linear Region

Looking at Figure 20, the purpose of this section is to make sure that ^{−1} on the locus of points, expressed by Eq. 35, and indicated via a dashed brown curve with magenta diamonds in plot (a), that such

Enforcing the inequality Eq. 55, and assuming a positive polarity for the offset current according to

the function

the characteristic

establishes the requirement for the point of intersection between the

Given that, with inequality Eq. 55 holding true, Eq. 58 expresses the condition under which ^{−1} in the linear region, the phase plane trajectories, lying, therein, below (above) the

##### 5.2.2 Recall M-CNN Cell DRM_{2} Synthesis in the Saturation Regions

Referring to Figure 20, the intention of this section is to establish the existence of a ^{−1} locus in the negative (positive) saturation region, as expressed by Eqs. 34, 36, and highlighted via a dashed brown curve with magenta diamonds in plot (b) ((c)), and to ensure that it forms, together with the ^{−1} locus, defined in Eqs. 31, 33, a stable equilibrium point, as given in Eqs. 39, 45, and shown as a black filled circle in plot (b) ((c)).

As anticipated earlier, mathematical calculations reveal that the

Focusing now on the domain

ensures that the graph of the function ^{−1} locus, the equilibrium Eq. 45, shown as a black circle on the same plot. On the basis of the behaviour of the vector field in the negative (positive) saturation region, as established by conditions Eqs. 48, 49 and Eqs. 50, 52, phase-plane trajectories below or above the ^{−1} locus Eqs. 34, 36 are bound to move north-eastward (north-westward) or south-eastward (south-westward), visiting the yellow II (cyan III) or gray IV (green I) regions, as qualitatively sketched in the viewgraph of Figures 20B,C, unvealing the local stability of the equilibrium point Eqs. 39, 45, as indicated through the filling structure of its black circle symbol.

All in all, the application of our stepwise system-theoretic M-CNN design method to the cell model Eqs. 18, 19, identifies five inequalities, specifically Eqs. 55, 65–68. Replacing the linear resistor in parallel to the capacitor in the memcomputing core of the cell circuit of Figure 8 with an open circuit^{43}, the system of five inequalities may be reduced to the triplet of non-redundant conditions Eqs. 55, 66, and 67. Solving them through the geometric analysis method, all the points lying in the magenta region of the *z*–_{2}-centered M-CNN design methodology. The binary image, shown in plot (a), illustrates graphically the resistive states of all the memristors. A white (black) pixel in a given position of this image reveals that the memristor in the equivalent location of the RECALL M-CNN stores the lowest (highest) possible resistance before the memory reading task is initiated. The image in plot (b) visualizes through a uniform gray color of appropriate tone^{44} the common initial condition assigned to each capacitor voltage, set to

**FIGURE 21**. **(A)** Black-and-white checkerboard visualizing the binary data stored in the locally distributed memristive memory bank previous to their retrieval. A white (black) pixel at row *i* and column *j* in this image reveals the on (off) resistance of the memristor in the M-CNN cell **(B)** Grey-scale image indicating the common **(C)** Output binary image coding the steady-state output voltages of all the RECALL M-CNN processing elements. A black (white) pixel in correspondence to the

## 6 Discussion

The theory presented in this work is independent of the memristor model adopted in M-CNN circuit design. In fact, it is a general theory, which may be applied to a much wider range of nonlinear dynamical circuits other than the cellular arrays analyzed in the manuscript. It is worth to pinpoint that the Second-Order Dynamic Route Map (DRM_{2}), around which the design methodology proposed in this paper is centered, extends the classical Dynamic Route Map (DRM) (Chua, 1998; Chua and Roska, 2002), applicable to first-order systems only, allowing to draw a complete picture of the local and global behaviour of any second-order dynamical system. For example, it shall constitute the system-theoretic tool of reference for a thorough study of the nonlinear dynamics of memristive devices with two state variables*.* In this work the DRM_{2} is adopted to investigate the spatio-temporal phenomena emerging in each of the second-order memristive cells of a two-dimensional cellular array. The two degrees of freedom of each cell are the voltage across a linear capacitor and the state of a first-order non-volatile memristor. The model (Pershin et al., 2009) adopted for the resistance switching memory in this work is a simple yet accurate mathematical description of a physical memristor realization (Jo et al., 2009). The reason behind the choice of this particular piecewise-linear memristor model for our study lies behind the pedagogical nature of this paper. In other words, since the aim of this manuscript is to provide researchers with powerful system-theoretic methods to analyze memristive cellular arrays, we found it useful to adopt an analytically tractable memristor model, in order to allow the determination of closed-form analytical expressions for the nullclines as well as for the equilibria of each State Dynamic Portrait (SDP) of a given DRM_{2}. The systematic M-CNN design methodology, presented in this paper, allows to derive optimal values for the cell parameters of each second-order cell _{2} would undergo inevitable changes, since the shape of the nullclines, the number, position, and stability of the second-order cell equilibria in the memristor state-capacitor voltage phase plane, the rules dictating how the sign of the time derivatives of the two states change across the

## 7 Conclusion

The motto “linearize-then-analyze”, which electrical engineers have been advocating for generations, should not drive the investigation of highly nonlinear memristive devices, circuits and systems, which are being developed in our times through disruptive nanotechnologies with the intention to foster progress in integrated circuit (IC) design beyond the Moore era. In fact, given that linear analysis techniques are unable to gain a deep insight into the behaviour of a nonlinear system, the availability of a partial picture of the dynamics of a novel nano-device prevents its conscious use in IC design. Recurring to nonlinear system theory is thus absolutely necessary to unfold the full potential of memristors in electronics. However, the conversion of classical circuits to memristive equivalents might require the adaptation of classical nonlinear system-theoretic analysis and design techniques, as is the case in this study. Cellular Nonlinear Networks (CNNs) (Chua and Yang, 1988a; Chua and Yang, 1988b) constitute one of the earliest examples of a non-von Neumann computing architecture, where data processing and storage tasks are locally distributed across a multi-dimensional array of locally coupled dynamical systems. In analogue hardware implementations of these bio-inspired computing structures, the cells typically feature one degree of freedom. As a result, the Dynamic Route Map (DRM) graphical tool, a powerful system-theoretic technique for the analysis of first-order systems, is applicable to gain a full understanding of the dynamics of these cells. Further, a rigorous procedure, employing the DRM analysis method, and leading to the derivation of an optimal solution for an inequality set (IS), which constrain number, and stability of cell equilibria for each of the possible combinations of inputs and/or initial conditions, allows to program the cellular network for the robust execution of a predefined computing task. The adoption of memristors in new designs of cell and coupling circuitry may allow to extend the processing functionalities and/or the computing efficiency of traditional dynamic arrays, thanks to the enrichment of the spectrum of dynamical phenomena, which may emerge within the cellular medium, while allowing to improve the spatial resolution of CNN analogue hardware realisations, concurrently. It is thus timely to investigate the impact of the introduction of memristors in new CNN designs. This work consider a first class of Memristor CNNs (M-CNNs), in which a first-order non-volatile resistance switching memory is inserted in parallel to the capacitor in each cell of a two-dimensional time- and space-invariant standard CNN. Given that the cells of each M-CNN from the proposed class feature two degrees of freedom, the DRM analysis methodology is no longer pertinent to gain insight into their data processing capabilities. A novel graphical tool, inspired to the Phase Portrait concept (Strogatz, 2000) from the theory of nonlinear dynamics, constituting the natural extension of the classical DRM system-theoretic technique to dynamical systems with two degrees of freedom, and called Second-Order Dynamic Route Map (DRM_{2}) (Tetzlaff et al., 2020), may allow to gain a deep insight into the dynamical phenomena emerging in cellular arrays with second-order memristive cells (Ascoli et al., 2020b), enabling to draw, finally, a codimension-2 bifurcation diagram, referred to as M-CNN Primary Mosaic, which specifies all the possible stable and unstable equilibria, which a cell may admit for each combination of self-feedback synaptic weight _{2} graphical tool, and leading to the derivation of an optimal solution of an IS, which shape the phase portrait of each cell in such a way that solutions of the CNN model equations may approach predefined equilibria for each of the possible combinations of inputs and initial conditions, allows to tune the parameters of the cellular array for a variability-tolerant accomplishment of a prescribed signal processing task or of a predefined memory operation. This work, contributing to the establishment of solid foundations of M-CNN theory, highlights the huge potential of memristive mem-processing structures for edge computing applications, and is expected to serve as a source of inspiration for future studies intended to verify the theoretical predictions on the beneficial impact of resistance switching memories on the performance of cellular nonlinear arrays.

## Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

## Author Contributions

AA and RT conceived the main idea of the paper. AA developed the systematic M-CNN design methodology, made all the analytical calculations, run the complete set of simulations to confirm the theoretical derivations, and wrote the whole manuscript. RT, SK, and LC supported the research with inspiring suggestions, precious guidelines, and insightful advices.

## Funding

This work has been partially supported by the Czech Science Foundation under grant No. 18–21608S, Czech Republic. LC is supported in part by AFOSR Grant No. FA 9550-18-1-0016.

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Footnotes

^{1}*In-memory computing* refers to the partial/temporary use of data storage units for information processing purposes, while *mem-computing* is associated to the adoption of computing systems to perform memory read/write operations on demand, as is the case for the processing elements of the proposed memristive cellular array.

^{2}Interestingly, it has been recently shown (Zhang et al., 2020) that the Cardiac Purkinje Fiber (CPF), which is the last branch of the heart conduction system, may be described via a modified variant of the memristive Hodgkin-Huxley equations, revealing the ubiquituous presence of memristors in living cells.

^{3}A CNN cell, exhibiting physical couplings to the eight closest processing elements, is said to have a sphere of influence of unitary radius, or, alternatively, a 3×3 local neighborhood.

^{4}The determination of a numerical solution for the *M*×*N* ODEs, governing the time evolution of the states of all the processing elements, requires the preliminary assignment of an *initial condition**boundary conditions* (Chua and Roska, 2002), fixing the input voltage *virtual cell*. With reference to a two-dimensional CNN, in which each cell exhibits a sphere of influence of unitary radius, a processing element

^{5}It is important to note that only CNN hardware realizations with such a basic coupling configuration have been developed so far (Vázquez et al., 2018).

^{6}Since the CNN is space-invariant, from now onwards the following assumptions are made:

^{7}As reported in the template library (Karacs et al., 2018), a very large class of computing tasks may be executed on the basis of the stable equilibria, the CNN states dynamically evolve to, or, more precisely, of the respective cell output voltages.

^{8}Images to be processed, consisting of *M* rows and *N* columns. The output voltages, that the cells exhibit at equilibrium, and that constitute the result of a given computing task, are mapped onto an output image for visualization purposes. Regarding the correspondence between the color of the image pixel at row *i* and column *j* and the real value of the respective CNN cell input voltage

^{9}The viewgraphs in Figure 3 have been derived under the following parameter setting:

^{10}In the case

^{11}The effect of the offset current is to shift the Internal DP Characteristic, explaining the origin for the name of the

^{12}This theorem holds also for standard space-variant CNNs.

^{13}For

^{14}A standard space-invariant CNN, in which the cell *uncoupled* (Chua and Roska, 2002) if each feedback synaptic weight *non-autonomous*. Each rule in a suitable set, which an uncoupled nonlinear dynamic array is obliged to obey, so as to execute a preliminarily specified data processing task, dictates the equilibrium *local*. In case the aforementioned CNN is coupled, on the other hand, some of the rules – referred to as global – for the cell

^{15}In the theory of nonlinear dynamics, a quantitative change in the behaviour of a system, occurring during the sweep of a control parameter, is referred to as a *bifurcation* phenomenon (Strogatz, 2000).

^{16}It is instructive to observe that a very large number of fundamental image processing operations are possible adopting the class of standard space-invariant uncoupled CNNs, as may be inferred by inspecting the template library (Karacs et al., 2018). For each CNN, belonging to this class and processing still images, the cell offset current from Eq. 13 is always a constant, and the respective Shifted DP Characteristic is invariant over time. As anticipated earlier, with reference to a coupled CNN, which process still images, a cell, which satisfies the bistability condition Eq. 12, typically features a Shifted DP Characteristic, which continually moves vertically up or down until the time instant

^{17}The values of some of the parameters in the cell circuit of Figure 2 are fixed as reported in the caption of Figure 5B. Since the bistability condition Eq. 12 holds true, the slope of the piecewise-linear locus in the linear region–refer to Eq. 9 is strictly positive.

^{18}In case *semistable equilibrium*, which attracts only trajectories, which are initiated from one of its two sides.

^{19}In the worst-case scenario from rule 3 the cell is found to be bistable provided

^{20}In case

^{21}Typically, the Full-Range model (Vázquez et al., 1993) is used in place for the Chua-Yang mathematical description of Eq. (1) to limit the range of admissible values for the cells’ states, thus simplifying the hardware realisation of the CNN paradigm.

^{22}In state-of-the-art sensor-processor arrays the input to the processing element, lying in correspondence to the

^{23}In the class of M-CNNs, investigated in this thesis, memristors are employed only for the design of the cell circuit. Their use for the circuit implementation of the synaptic couplings shall be the focus of future studies.

^{24}Throughout this chapter the state, voltage, and current of the memristor

^{25}The numerical integration of the *initial condition*

^{26}The Pershin and Di Ventra model from (Pershin et al., 2009) adopts the Heaviside functions

^{27}A phase plane trajectory is the locus of points

^{28}This 19-parameter set is often referred to as *gene* in CNN theory (Chua and Roska, 2002), since its 19 elements crucially affect the spatio-temporal phenomena, which may emerge in the space-invariant array of locally coupled cells, similarly as the DNA genetic content has a significant impact on the dynamical evolution of living beings. The 9 feedback (feedforward) synaptic weights in the set **A** (**B**), referred to as *feedback (feedforward) template* in CNN theory (Chua and Roska, 2002).

^{29}The conditions, under which a space-invariant M-CNN is uncoupled, are the same as defined in section 2.3 for a space-invariant CNN: *local* (Chua and Roska, 2002), i.e., in general, they depend upon the input voltage *global* (Chua and Roska, 2002)–for the cell

^{30}The M-CNN computing paradigm, this section is focused upon, revolves around the asymptotic convergence of the state vector

^{31}Despite template optimization (Chua and Roska, 2002) does not constitute the focus of this research study, in some cases, in order to improve the robustness of the M-CNN design, it may be adviceable to include additional inequalities, which may endow the design with a good tolerance to parameter variability. As examples, one could enforce a minimal distance between certain

^{32}These very same values will be assigned to all the invariable cell circuit parameters in each of the M-CNN designs discussed below.

^{33}Indicating the memristor resistance

^{34}As will be clarified shortly, the scenario

^{35}Given that, as discussed in section 3.4, under

^{36}The direction of motion of the state vector

^{37}The function of the linear resistor is to decrease the absolute value of the

^{38}As anticipated earlier, for each combination of indices

^{39}Despite, theoretically, under the hypothesis of either rule, the two-dimensional state vector *t* tends to

^{40}A standard space-invariant uncoupled CNN, in which each cell features a *isolated* (Chua and Roska, 2002) if each feedforward synaptic weight

^{41}Similarly as in the EDGE M-CNN design, the addition of a linear resistor in parallel to the capacitor within the memcomputing core of the processing element circuit of Figure 8 allows to keep within reasonable limits the modulus of the voltage, falling across the resistance switching memory at equilibrium, for each of the two possible cell input voltage values. This preventive measure is of particular importance in view of a future hardware realization of the bio-inspired memristive array under study.

^{42}The task of the RECALL M-CNN may be considered accomplished as soon as the outputs of all its processing elements attain their final values. This occurs at the steady-state time instant

^{43}The modulus of the voltage, falling across the cell memristor at equilibrium, was found to be reasonably small, even for

^{44}The graphical illustration convention, adopted here for visualizing each initial cell capacitor voltage is analogous as the one established in the discussion of the EDGE M-CNN design: the closer is its real value, lying in the set

## References

Ascoli, A., Tetzlaff, R., Kang, S. M, and Chua, L. O. (2020a). Theoretical foundations of memristor cellular nonlinear networks: a DRM-based method to design memcomputers with dynamic memristors. *IEEE Trans. Circuits Systems–I: Regular Pap.* 67 (8), 2753–2766. doi:10.1109/tcsi.2020.2978460

Ascoli, A., Messaris, L., Tetzlaff, R., and Chua, L. O. (2020b). Theoretical foundations of memristor cellular nonlinear networks: stability analysis with dynamic memristors. *IEEE Trans. Circuits Systems–I: Regular Pap.* 67 (4), 1389–1401. doi:10.1109/tcsi.2019.2957813

Biolek, Z., Biolek, D., and Biolkova, V. (2009). Spice model of memristor with nonlinear dopant drift. *Radioengineering* 18 (2), 210–214.

Bohaichuk, S. M., Kumar, S., Pitner, G., McClellan, C. J., Jeong, J., Samant, M. G., et al. (2019). Fast spiking of a mott VO – carbon nanotube composite device. *Nano Lett.* 19, 6751–6755. doi:10.1021/acs.nanolett.9b01554

L. O. Chua (Editor) (1998). *CNN: a paradigm for complexity. World scientific series on nonlinear science, series A, vol. 31*. Singapore: World Scientific.

Chua, L. O. (2015). Everything you wish to know about memristors but are afraid to ask. *Radioengineering* 24 (2), 319–368. doi:10.13164/re.2015.0319

Chua, L. O. (2018a). Five non-volatile memristor enigmas solved. *Appl. Phys. A*. 124 (8), 563. doi:10.1007/s00339-018-1971-0

Chua, L. O., Sbitnev, V., and Kim, K. (2012). Hodgkin-Huxley axon is made of memristors. *Int. J. Bifurcation Chaos* 22 (3), 1230011. doi:10.1142/s021812741230011X

Chua, L. (2014). If It’s Pinched, It’s a Memristor. *Semicond. Sci. Technol. Spec. Issue. Memristive Devices* 29 (10), 42. doi:10.1088/0268-1242/29/10/104001

Chua, L. O., and Kang, S. M. (1976). Memristive devices and systems. *Proc. IEEE* 64 (2), 209–223. doi:10.1109/proc.1976.10092

Chua, L. O. (2005). Local activity is the origin of complexity. *Int. J. Bifurcation Chaos* 15 (11), 3435–3456. doi:10.1142/s0218127405014337

Chua, L. O. (1971). Memristor: the missing circuit element. *IEEE Trans. Circuit Theor.* 18 (5), 507–519. doi:10.1109/tct.1971.1083337

Chua, L. O. (2018b). Memristors: remembrance of things past. *IEEE Micro* 38 (5), 7–12. doi:10.1109/mm.2018.053631136

Chua, L. O. (2011). Resistance switching memories are memristors. *Appl. Phys. A.* 102, 765–783. doi:10.1007/s00339-011-6264-9

Chua, L. O., and Roska, T. (2002). *Cellular neural networks and visual computing: foundations and applications*. 1 edition. Cambridge, England: Cambridge University Press.

Chua, L. O., and Yang, L. (1988b). Cellular neural networks: applications. *IEEE Trans. Circuits Systems–i* 35 (10), 1273–1290. doi:10.1109/31.7601

Chua, L. O., and Yang, L. (1988a). Cellular neural networks: theory. *IEEE Trans. Circuits Systems–I* 35 (10), 1257–1272. doi:10.1109/31.7600

Di Marco, M., Forti, M., and Pancioni, L. (2017a). Memristor standard cellular neural networks computing in the flux-charge domain. *Neural Netw.* 93, 152–164. doi:10.1016/j.neunet.2017.05.009

Di Marco, M., Forti, M., and Pancioni, L. (2017b). Convergence and Multistability of Nonsymmetric Cellular Neural Networks With Memristors. *IEEE Trans. Cybern.* 47 (10), 2970–2983. doi:10.1109/TCYB.2016.2586115

Di Marco, M., Forti, M., and Pancioni, L. (2018). New Conditions for Global Asymptotic Stability of Memristor Neural Networks. *IEEE Trans. Neural Netw. Learn. Syst.* 29, 1822–1834. doi:10.1109/TNNLS.2017.2688404

Duan, S., Hu, X., Dong, Z., Wang, L., and Mazumder, P. (2015). Memristor-Based Cellular Nonlinear/Neural Network: Design, Analysis, and Applications. *IEEE Trans. Neural Netw. Learn. Syst.* 26 (6), 1202–1213. doi:10.1109/TNNLS.2014.2334701

Fujitsu Ltd (2019). Fujitsu announcement. Available at: https://www.fujitsu.com/global/products/devices/semiconductor/memory/reram/ (Accessed July 30, 2019).

GlobalFoundries Ltd (2018). GlobalFoundries announcement. Available at: https://spectrum.ieee.org/nanoclast/semiconductors/devices/globalfoundries-halts-7nm-chip-development (Accessed August 28^{th}, 2018).

Hu, M., Strachen, J. P., Li, Z., Grafals, E. M., Davila, N., Lam, S., et al. (2016). “Dot-product engine for neuromorphic computing: programming 1T1M crossbar to accelerate vector-matrix multiplication,” in Proceedings of the ACM/EDAC/IEEE 53rd Annual Design Automation Conference (DAC), Austin, TX, USA, June 5–9, 2016.

Itoh, M., and Chua, L. O. (2003). Designing CNN genes. *Int. J. Bifurcations Chaos* 13 (10), 2739–2824. doi:10.1142/s0218127403008375

Jo, S. H., Kim, K. H., and Lu, W. (2009). High-density crossbar arrays based on a Si memristive system. *Nano Lett.* 9 (2), 870–874. doi:10.1021/nl8037689

Karacs, K., Cserey, Gy., Zarándy, Á., Szolgay, P., Rekeczky, Cs., Kék, L., et al. (2018). CNN template library. Available at: http://cnn-technology.itk.ppke.hu/Template_libraryv4.0beta2.pdf (Accessed April 19, 2021).

Kumar, S., Williams, R. S., and Wang, Z. (2020). Third-order nano-circuit elements for neuromorphic engineering. *Nature* 585, 518–523. doi:10.1038/s41586-020-2735-5

Landauer, R. (1988). Dissipation and noise immunity in computation and communication. *Nature* 335 (6193), 779–784. doi:10.1038/335779a0

Li, C., Hu, M., Li, Y., Jiang, N., Ge, N., Montgomery, E., et al. (2018). Analogue signal and image processing with large memristor crossbars. *Nat. Electron.* 1, 52–59. doi:10.1038/s41928-017-0002-z

Moore, G. E. (1965). Cramming more components onto integrated circuits. *Electronics* 38 (8), 114–117. doi:10.2307/3756714

Panasonic Ltd. (2013). Panasonic announcement. Available at: https://news.panasonic.com/global/press/data/2013/07/en130730-2/en130730-2.html (Accessed July 30^{th}, 2013).

Pershin, Y. V., and Di Ventra, M. (2011). Memory effects in complex materials and nanoscale systems. *Adv. Phys.* 60 (2), 145–227. doi:10.1080/00018732.2010.544961

Pershin, Y. V., Fontaine, S. L., and Di Ventra, M. (2009). Memristive model of amoeba learning. *Phys. Rev. E*. 89, 021926. doi:10.1103/PhysRevE.80.021926

Roska, T. (1993). The CNN universal machine: an analogic array computer. *IEEE Trans. Circuits Syst.–II: Analog Digital Signal Process.* 40 (3), 163–173. doi:10.1109/82.222815

Strogatz, S. H. (2000). *Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering (studies in nonlinearity)*. 1 edition. Florida, US: CRC Press.

Tetzlaff, R., Ascoli, A., Messaris, I., and Chua, L. O. (2020). Theoretical foundations of memristor cellular nonlinear networks: memcomputing with bistable-like memristors. *IEEE Trans. Circuits Syst.–I: Regular Pap.* 67 (2), 502–515. doi:10.1109/tcsi.2019.2940909

Toshiba Ltd. (2012). Smart photo sensor SPS02. Available at: https://www.toshiba-teli.co.jp/pdf/industrial-camera-en/sps02e.pdf (Accessed February 1st, 2012).

Vázquez, A. R., Espejo, S., Domínguez-Castro, R., Huertas, J. L., and Sánchez-Sinencio, E. (1993). Current-mode techniques for the implementation of continuous- and discrete-time cellular neural networks. *IEEE Trans. Circuits Syst. II: Analog Digit. Signal Process.* 40 (3), 132–146. doi:10.1109/82.222812

Vázquez, A. R., Fernández-Berni, J., Leñero-Bardallo, J. A., Vornicu, I., and Carmona-Galán, R. (2018). CMOS vision sensors: embedding computer vision at imaging front-ends. *IEEE Circuits Syst. Mag.* 18 (2), 90–107. doi:10.1109/MCAS.2018.2821772

Weiher, M., Herzig, M., Tetzlaff, R., Ascoli, A., Mikolajick, T., and Slesazeck, S. (2019). Pattern formation with local active S-type NbOx memristors. *IEEE Trans. Circuits Syst.–I: Regular Pap.* 66 (7), 2627–2638. doi:10.1109/tcsi.2019.2894218

Williams, R. S. (2017). What’s next? [The end of Moore’s law]. *IEEE Comput. Sci. Eng.* 19 (2), 7–13. doi:10.1109/mcse.2017.31

Zarándy, Á. (2003). The art of CNN template design. *Int. J. Circuit Theor. Appl.* 27 (1), 5–23. doi:10.1002/(sici)1097-007x(199901/02)27:1<5::aid-cta38>3.0.co;2-c

Keywords: memristor, bio-inspired mem-computing machines, cellular nonlinear networks, nonlinear circuit theory, nonlinear system theory

Citation: Ascoli A, Tetzlaff R, Kang S-MS and Chua L (2021) System-Theoretic Methods for Designing Bio-Inspired Mem-Computing Memristor Cellular Nonlinear Networks. *Front. Nanotechnol.* 3:633026. doi: 10.3389/fnano.2021.633026

Received: 24 November 2020; Accepted: 25 January 2021;

Published: 12 May 2021.

Edited by:

Huanglong Li, Tsinghua University, ChinaReviewed by:

Carlos Sánchez-López, Autonomous University of Tlaxcala, MexicoXiaoping Wang, Huazhong University of Science and Technology, China

Copyright © 2021 Ascoli, Tetzlaff, Kang and Chua. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Alon Ascoli, alon.ascoli@tu-dresden.de