Memristor Circuits for Simulating Neuron Spiking and Burst Phenomena

Since the introduction of memristors, it has been widely recognized that they can be successfully employed as synapses in neuromorphic circuits. This paper focuses on showing that memristor circuits can be also used for mimicking some features of the dynamics exhibited by neurons in response to an external stimulus. The proposed approach relies on exploiting multistability of memristor circuits, i.e., the coexistence of infinitely many attractors, and employing a suitable pulse-programmed input for switching among the different attractors. Specifically, it is first shown that a circuit composed of a resistor, an inductor, a capacitor and an ideal charge-controlled memristor displays infinitely many stable equilibrium points and limit cycles, each one pertaining to a planar invariant manifold. Moreover, each limit cycle is approximated via a first-order periodic approximation analytically obtained via the Describing Function (DF) method, a well-known technique in the Harmonic Balance (HB) context. Then, it is shown that the memristor charge is capable to mimic some simplified models of the neuron response when an external independent pulse-programmed current source is introduced in the circuit. The memristor charge behavior is generated via the concatenation of convergent and oscillatory behaviors which are obtained by switching between equilibrium points and limit cycles via a properly designed pulse timing of the current source. The design procedure takes also into account some relationships between the pulse features and the circuit parameters which are derived exploiting the analytic approximation of the limit cycles obtained via the DF method.

1. Introduction

The ever-growing need of computing power to handle data intensive applications is seriously challenging conventional digital von Neumann computing architectures (Bonomi et al., 2012; Satyanarayanan, 2017; Williams, 2017). The physical separation between the computing and memory units can indeed generate long latency time and large energy consumption when data intensive tasks are performed. In this context, researchers look at the emerging nanoscale analog devices, such as memristors, as a viable approach for developing new computing paradigms, based on in-memory and analog computation, which are potentially capable to overcome the limitations of conventional computer architectures (Waldrop, 2016; Zidan et al., 2018; Krestinskaya et al., 2020).

The memristor (a shorthand for memory resistor) is the fourth basic passive circuit element theoretically introduced by Prof. Leon Chua in 1971 (Chua, 1971). The appealing features of a memristor are non-volatility, i.e., the memristor capability to hold data in memory without the need of a power supply, and non-linearity, which makes memristor circuits capable to generate quite a rich variety of oscillatory and complex dynamics. The combination of these two features enables in-memory computing, i.e., the co-location of computation and memory in the same device (Ielmini and Wong, 2018). In-memory computing, which resembles a basic principle of brain computation, can provide several advantages, such as low energy consumption, high bandwidths, excellent scalability, thus lending itself as quite a promising novel computing approach in the field of artificial intelligence and big data (Ielmini and Pedretti, 2020). Memristor circuits have been already positively used to address several tasks, including the solution of global optimization, constraint satisfaction and linear algebra problems (Yang et al., 2013; Wang et al., 2015; Kumar et al., 2017). They are also used as building blocks in reservoir computing systems (Du et al., 2017) and neuromorphic computing for on-line signal processing tasks (Di Marco et al., 2016, 2017; Di Marco et al., 2017; Ascoli et al., 2020a,b).

Since the very beginning it was clear that understanding the peculiar dynamical features of memristor circuits is the key step for developing analog in-memory computing schemes. It has been definitely shown that memristor circuits are capable to generate quite a large variety of dynamical behaviors (Corinto et al., 2011; Corinto et al., 2019; Pershin and Di Ventra, 2011; Xu et al., 2016; Yuan et al., 2016; Di Marco et al., 2018), including bursting oscillations (see e.g., Wang et al., 2019 and references therein). Recently, a new technique, named flux-charge analysis method (FCAM), has been introduced for the analysis of memristor circuits in the flux-charge domain (Corinto and Forti, 2016, 2017, 2018). FCAM permits to show that the dynamical richness displayed by memristor circuits is a consequence of the property that the state space of any given circuit, i.e., with its parameter having fixed values, contains infinitely many invariants manifolds (foliation property of the state space). This specific property implies that in memristor circuits there is the coexistence of infinitely many different attractors, a property referred to as multistability. Moreover, structural changes of the attractors are observed when the initial conditions are varied while keeping constant the values of the circuit parameters, a peculiar phenomenon which is referred to as bifurcations without parameters (Corinto and Forti, 2017; Di Marco et al., 2018; Innocenti et al., 2019b). In particular, Di Marco et al. (2018), Innocenti et al. (2019b) have employed techniques within the Harmonic Balance (HB) context for predicting limit cycles and their bifurcations by first showing that the dynamics of the memristor circuit admits an equivalent input-output representation, which has been recently extended also to circuits containing memory elements (Innocenti et al., 2020).

Among other applications, it has been soon recognized that memristors can be successfully employed as synapses in neuromorphic circuits (see e.g., Jo et al., 2010; Adhikari et al., 2012; Thomas, 2013; Kim et al., 2015; Hu et al., 2016; Hong et al., 2019). It has been also pointed out that memristor circuits appear to be suited for modeling some features of the dynamics of neurons. Some contributions provide a memristor representation of popular neuron models (see e.g., Chua et al., 2012; Usha and Subha, 2019 for the Hodgkin-Huxley axon), while others show how to mimic some typical dynamics displayed by cortical neurons (see e.g., Babacan et al., 2016; Nakada, 2019, and references therein). In Innocenti et al. (2019a), it is shown that dynamics of the memristor Murali-Lakshmanan-Chua circuit, equipped with an independent pulse programmed input source and simple comparator and hysteresis feedback blocks, can resemble some dynamical behaviors of cortical neurons. Finally, it is worth noting that HB techniques have been suitably employed for the analysis of neural oscillations (see e.g., Chen et al., 2008; Matsuoka, 2011 and references therein).

The purpose of this paper is to show how memristor circuits can be exploited for modeling some features of the neurons dynamics (see e.g., Izhikevich, 2000, 2007). The basic ideas are to exploit multistability, i.e., the fact that a memristor circuit contains infinitely many attractors (equilibria, limit cycles, …), and to employ a suitable pulse-programmed input for controlling multistability, i.e., for switching among the different attractors. Controlling multistability is currently a research field of general interest (see e.g., Pisarchik and Feudel, 2014 and references therein) and some contributions to the problem of steering the memristor circuit dynamics toward the attractors contained in one of the invariant manifolds have been given (Chen et al., 2018; Corinto and Forti, 2018; Corinto et al., 2019; Di Marco et al., 2019, 2020a,b). In particular, Di Marco et al. (2020b,a), Di Marco et al. (2021) have shown that the dynamics of a second order R, L, C circuit connected with a charge-controlled memristor can be steered, via a pulse programmed source, toward a pre-assigned invariant manifold where it converges toward the attractor contained in the manifold itself. In this paper, such a simple memristor circuit is used to mimic some aspects of the dynamics displayed by cortical neurons in response to an external stimulus. Section 2 describes the circuit and formulates the problem of interest. Specifically, the memristor charge should exhibit some typical neuron dynamics when the current source is suitably pulse-programmed. It is assumed that the pulses and their timing are generated by some given hardware mechanisms, whose design is not the object of the paper. Section 2 also illustrates the foliation property of the circuit state space in the input-less case, by characterizing the infinitely many invariant manifolds and the differential equations governing the dynamics onto them. Section 3 investigates the dynamics of the input-less circuit by showing that each manifold contains as attractor either a stable equilibrium point or a stable limit cycle. The limit cycles analysis is performed via the Describing Function (DF) method, a classical technique within the HB context. By exploiting the coexistence of stable equilibrium points and limit cycles and the knowledge of the first order periodic approximations, section 4 provides a procedure for the design of a pulse timing of the current source ensuring that the memristor charge mimics some dynamical features of the neuron response. The sought behavior of the memristor charge is obtained via a suitable concatenation of convergent and oscillatory behaviors, which are generated by switching between different attractors according to the designed pulse timing. The relation between the pulse timing parameters and the circuit parameters is also discussed. Some conclusions are finally drawn in section 5.

2. Problem Formulation and Preliminaries

The aim of the paper is to show that the coexistence of many different attractors in memristor circuits permits to mimic some features of the dynamics of cortical neurons. Specifically, we consider the simple circuit depicted in Figure 1 which contains a resistor R, an inductor L, a capacitor C and a non-linear charge-controlled memristor.

FIGURE 1

Figure 1. Charge-controlled memristor circuit with R, L,C components and independent current source I_s.

The capacitor voltage, the inductor current, the memristor voltage and current are denoted by v_C, i_L, v_M, and i_M, respectively. The circuit is subject to an independent current source I_s, which is the input of the circuit. The charge-flux characteristic $\widehat{\phi }:ℝ\to ℝ$ of the memristor relating the charge q_M and the flux φ_M is assumed to have both a linear and a cubic term. Specifically,

$\begin{array}{ll}{\phi }_M=\widehat{\phi }(q_M)=-s_0{q}_M+\frac{s_1}{3}{q}_M^3,& (1)\end{array}$

where the constant terms s₀ and s₁ satisfy

$\begin{array}{ll}{s}_0>R,s_1>0.& (2)\end{array}$

Throughout the paper, we assume that all the circuit parameters R, L, C, s₀, s₁ have normalized values. Also, we consider arbitrary units for the time.

We want to show that the circuit is able to generate some characteristic dynamical behaviors of a neuron in response to a pulse stimulus, as the typical one depicted in Figure 2A. Specifically, the memristor charge q_M should display the time behavior of Figure 2B when the input source I_s provides a suitable stimulus.

FIGURE 2

Figure 2. (A) Schematic of an ideal voltage-pulse of the cellular membrane (action potential). (B) Reference shape of the memristor charge-pulse considered in this paper. For the sake of simplicity, the depolarization is assumed to occur without distinction between the sub- and the super-threshold branches, and the hyperpolarization dynamics admits small ripples during the convergence to the resting state.

It can be readily verified that the circuit dynamics is described by the following third-order system of differential equations

$\begin{array}{ll}\{\begin{array}{c}{D}{v}_C(t)=\frac{1}{C}\left(I_s(t)-i_L\right)\\ {D}{q}_M(t)=i_L\\ {D}{i}_L(t)=\frac{1}{L}\left(v_C+(s_0-R)i_L-s_1{q}_M^2{i}_L\right)\end{array}& (3)\end{array}$

where ${D}$ is the time-derivative operator¹ and v_C, q_M, and i_L are the state variables.

Let us consider the input-less case, i.e., I_s = 0. It has been shown (see e.g., Corinto and Forti, 2017) that memristor circuits enjoy the so-called foliation property, i.e., the memristor state space is decomposed into infinitely many invariant manifolds. The verification of this property for the considered circuit is readily obtained. Since I_s = 0, the first equation of Equation (3) can be rewritten equivalently as

$\begin{array}{ll}{D}{v}_C(t)+\frac{1}{C}{i}_L(t)={D}{v}_C(t)+\frac{1}{C}{D}{q}_M(t)=\frac{1}{C}{D}(q_M(t)+C{v}_C(t))=0& (4)\end{array}$

which means that the quantity q_M(t) + Cv_C(t) is constant along the solutions of Equation (3) with initial conditions v_C(t₀), q_M(t₀), and i_L(t₀) at time t = t₀, i.e.,

$\begin{array}{ll}{q}_M(t)+C{v}_C(t)=q_M(t_0)+C{v}_C(t_0)\forall t\ge t_0.& (5)\end{array}$

Hence, in the input-less case the state space of the memristor circuit is decomposed into infinitely many invariant manifolds of the form

$\begin{array}{ll}M(Q_0)=\{(v_C,q_M,i_L):q_M(t)+C{v}_C(t)=Q_0,\forall t\ge t_0\},& (6)\end{array}$

where Q₀ is the index of the manifold whose value depends on the circuit initial conditions according to the relation Q₀ = q_M(t₀) + Cv_C(t₀). Note that the invariant manifolds M(Q₀) are planar surfaces parameterized by the manifold index Q₀, as illustrated in Figure 3.

FIGURE 3

Figure 3. Invariant manifolds of the input-less circuit: the state space trajectory generated by the solution of Equation (3) with initial conditions v_C(t₀), q_M(t₀), i_L(t₀) (marked with ◦) at time t = t₀ belongs to the manifold M(Q₀) with Q₀ = q_M(t₀) + Cv_C(t₀) for all t ≥ t₀.

We have then established that in the input-less case the circuit dynamics is confined to lie onto one of the invariant manifolds. The differential equations governing the dynamics onto the invariant manifold M(Q₀) can be readily singled out. Since the first equation of Equation (3) has been used to obtain M(Q₀), it follows that the dynamics is characterized by the two remaining equations once v_C(t) is replaced by $\frac{1}{C}(Q_0-q_M(t))$. Hence, the dynamics onto M(Q₀) obeys the second-order system of differential equations

$\begin{array}{ll}\{\begin{array}{c}{D}{q}_M(t)=i_L(t)\\ {D}{i}_L(t)=\frac{1}{LC}(Q_0-q_M(t))+\frac{s_0-R}{L}{i}_L(t)-\frac{s_1}{L}{q}_M^2(t)i_L(t).\end{array}& (7)\end{array}$

Clearly, the complete dynamical behaviors of the input-less circuit can be obtained by collecting all the dynamics confined to lie onto each invariant manifold M(Q₀). The dynamical analysis of Equation (7) for any value of the index manifold Q₀ will be pursued in section 3. Note that Q₀ can be seen as a parameter of Equation (7) which however depends on the circuit initial conditions and hence it has a different nature with respect to the circuit parameters R, L, C, s₀, and s₁.

3. Input-Less Memristor Circuit Dynamics

In this section, we summarize the properties of the dynamics onto M(Q₀) by investigating system (7). Specifically, the equilibrium points and their stability features are dealt with in section 3.1, while the presence of limit cycles is considered in section 3.2.

3.1. Equilibrium Points

It readily follows that each invariant manifold M(Q₀) has a unique equilibrium point at (q_M, i_L) = (Q₀, 0). Taking into account (6), this implies that the circuit system (7) has an equilibrium point at (v_C, q_M, i_L) = (0, Q₀, 0) for any value of the index Q₀. Hence, the q_M-axis is a line of equilibrium points of the circuit.

To assess the stability properties of the equilibrium point of M(Q₀), we compute the Jacobian J(Q₀) of (7) at (q_M, i_L) = (Q₀, 0), getting

$\begin{array}{ll}J(Q_0)=\left(\begin{array}{cc}0& 1\\ -\frac{1}{LC}& \frac{s_1}{L}({\stackrel{-}{Q}}_0^2-Q_0^2)\end{array}\right),& (8)\end{array}$

where

$\begin{array}{ll}{\stackrel{-}{Q}}_0=\sqrt{\frac{s_0-R}{s_1}}.& (9)\end{array}$

The eigenvalues of J(Q₀) have a negative real part if $\left|Q_0\right|>{\stackrel{-}{Q}}_0$ and a positive real part if $\left|Q_0\right|<{\stackrel{-}{Q}}_0$, while they are pure imaginary at $Q_0=\pm {\stackrel{-}{Q}}_0$. This implies that the equilibrium point at (q_M, i_L) = (Q₀, 0) is ensured to be (locally) asymptotically stable if $\left|Q_0\right|>{\stackrel{-}{Q}}_0$ and unstable if $\left|Q_0\right|<{\stackrel{-}{Q}}_0$.

Hence, the equilibrium point of each manifold M(Q₀) with $\left|Q_0\right|>{\stackrel{-}{Q}}_0$ is an attractor of the input-less memristor circuit. Figure 4 summarizes the dynamical behavior around the equilibrium point of M(Q₀).

FIGURE 4

Figure 4. Stable (green half lines) and unstable (red segment) equilibrium points: the trajectories starting on the planes with $Q_0<-{\stackrel{-}{Q}}_0$ and $Q_0>{\stackrel{-}{Q}}_0$ converge toward the corresponding equilibrium points; the trajectory starting on the plane with $Q_0\in (-{\stackrel{-}{Q}}_0,{\stackrel{-}{Q}}_0)$ converges toward the stable limit cycle (green).

3.2. Limit Cycles

To investigate the presence of limit cycles on M(Q₀) we resort to the DF method, a classical technique within the HB context. The DF method has been widely used for analyzing oscillatory behaviors in non-linear feedback control systems (see e.g., Gelb and Vander Velde, 1968; Atherton, 1975; Mees, 1981; Khalil, 2002). Within the HB setting, other techniques have been proposed to predict bifurcations of limit cycles and more complex dynamics (Genesio and Tesi, 1992; Piccardi, 1994; Tesi et al., 1996; Basso et al., 1997; Bonani and Gilli, 1999; Di Marco et al., 2003; Innocenti et al., 2010). The first step to apply the DF method is to show that system (7) admits the representation of Figure 5 whose dynamics is governed by the following input-output relation

$\begin{array}{ll}{q}_M(t)=-L({D})\widehat{\phi }(q_M(t))+L_0({D})Q_0,& (10)\end{array}$

where $L({D})$ and $L_0({D})$ are suitable rational functions.

FIGURE 5

Figure 5. Equivalent input-output representation of system (7).

It is worth observing that this representation has an internal feedback interconnection between the linear subsystem $L({D})$ and the non-linear subsystem described by the memristor characteristic, while $L_0({D})$ is a feed-forward linear block driven by a constant input.

It can be readily verified that Equation (7) can be rewritten equivalently as a unique second order differential equation involving only q_M. Since $i_L={D}{q}_M$ the second equation of Equation (7) becomes

$\begin{array}{l}{{D}}^2{q}_M(t)=\frac{1}{LC}(Q_0-q_M(t))-\frac{R}{L}{D}{q}_M(t)\\ -\frac{1}{L}(-s_0{D}{q}_M(t)+s_1{q}_M^2(t){D}{q}_M(t)).& (11)\end{array}$

Taking into account that $-s_0{D}{q}_M(t)+s_1{q}_M^2(t){D}{q}_M(t)={D}{\widehat{\phi }}_M(t)$, the above differential equation can be rearranged in the following equivalent input-output form

$\begin{array}{ll}{q}_M(t)=-\frac{\frac{1}{L}{D}}{{{D}}^2+\frac{R}{L}{D}+\frac{1}{LC}}\widehat{\phi }(q_M(t))+\frac{\frac{1}{LC}}{{{D}}^2+\frac{R}{L}{D}+\frac{1}{LC}}{Q}_0.& (12)\end{array}$

Hence, system (7) can be equivalently described via the input-output relation (10) with $L({D})$ and $L_0({D})$ given by

$\begin{array}{ll}L({D})=\frac{\frac{1}{L}{D}}{{{D}}^2+\frac{R}{L}{D}+\frac{1}{LC}},L_0({D})=\frac{\frac{1}{LC}}{{{D}}^2+\frac{R}{L}{D}+\frac{1}{LC}}.& (13)\end{array}$

It is worth observing that $L({D})$ is exactly the equivalent admittance of the linear part of the circuit with I_S(t) = 0, i.e., the one seen at the terminals T′ and T^′′ in Figure 1 where the memristor is connected.

The DF method looks for a first-order approximation of a periodic output q_M(t) of the system of Figure 5, i.e.,

$\begin{array}{ll}{q}_M(t)=A+B\text{ cos }\omega t,B>0,\omega >0.& (14)\end{array}$

The corresponding periodic output φ_M(t) of period 2π/ω of the non-linear subsystem can be written as

$\begin{array}{l}\hfill {\phi }_M(t)=\widehat{\phi }(A+B\text{ cos }\omega t)=\\ N_0(A,B)A+N_1(A,B)B\text{ cos }\omega t+{\phi }_h(t),& (15)\end{array}$

where N₀(A, B) is the bias gain and N₁(A, B) is first harmonic gain, which are known as describing function terms, while φ_h(t) contains the higher order harmonics, i.e.,

$\begin{array}{ll}{\phi }_h(t)={\displaystyle \sum _{h=2}^{+\infty }}{\gamma }_h\text{ cos }(h\omega t+{\theta }_h)& (16)\end{array}$

with γ_h and θ_h, h = 2, …, being suitable constants.

The final step of the DF method consists in first computing the periodic signal of period 2π/ω given by the sum of the outputs y(t) and y₀(t) of the two linear subsystems driven by −φ_M(t) and the constant signal Q₀, respectively, then in equating the constant and the first order harmonic terms of the obtained periodic signal with the corresponding terms of q_M(t), i.e., A and B. Taking into account that the constant and first order harmonic terms of φ_M(t) can be rewritten in an exponential form as $N_0(A,B)A{e}^{0t}$ and $N_1(A,B)B(e^{j\omega t}+e^{-j\omega t})/2$, respectively, the periodic signal given by the sum of y(t) and y₀(t) can be computed by exploiting superposition and the expression of the response of a linear system which has the same exponential form of the input². The so computed periodic signal has the following expression

$\begin{array}{l}\hfill L(0)N_0(A,B)A+L_0(0)Q_0\hfill \\ +N_1(A,B)B\left(\frac{L(j\omega )e^{j\omega t}}{2}+\frac{L(-j\omega )e^{-j\omega t}}{2}\right)+y_h(t)& (17)\end{array}$

where y_h(t) contains the higher order harmonics generated by −φ_h(t). Then, by equating the constant term in Equation (17) with that of q_M(t), we get the real equation

$\begin{array}{ll}A+L(0)N_0(A,B)A=L_0(0)Q_0,& (18)\end{array}$

while by equating the first harmonic terms and taking into account that B cos ωt = B(e^jωt + e^−jωt)/2 we get the complex equation

$\begin{array}{ll}1+L(j\omega )N_1(A,B)=0.& (19)\end{array}$

In the HB approach, any first-order periodic signal (14) with A, B, and ω solving both Equations (18) and (19) is called a Predicted Limit Cycle (PLC) of the system of Figure 5. Also, it is expected that there exists a true limit cycle close to a PLC if the system of Figure 5 satisfies some conditions (see Atherton, 1975; Mees, 1981; Khalil, 2002 and references therein). These conditions basically rely on so-called filtering hypothesis which means that the non-linear subsystem does not generate large higher-order harmonics (i.e., φ_h(t) is small) and the two linear subsystems are low-pass filters (i.e., the gains |L(jhω)|, h = 2, …, are smaller than |L(0)| and |L(jω)|). Hence, this hypothesis requires that y_h(t) must be much smaller than the constant and first order harmonic terms of Equation (17), according to some quantitative measure. For instance, in Di Marco et al. (2018) the distortion index is used to this purpose.

Now, to apply the DF method to the system under consideration, we observe that from Equation (13) we get L(0) = 0 and L₀(0) = 1, while from Equations (1) and (15) it can be verified that the gains N₀(A, B) and N₁(A, B) have the following expressions:

$\begin{array}{ll}{N}_0(A,B)=-s_0+\frac{s_1}{3}\left(A^2+\frac{3}{2}{B}^2\right),& (20)\end{array}$

$\begin{array}{ll}{N}_1(A,B)=-s_0+s_1\left(A^2+\frac{1}{4}{B}^2\right).& (21)\end{array}$

Then, the Equation (18) reduces to

$\begin{array}{ll}A=Q_0,& (22)\end{array}$

and the complex equation (19) boils down to

$\begin{array}{ll}1+\frac{\frac{1}{L}j\omega }{-{\omega }^2+\frac{R}{L}j\omega +\frac{1}{LC}}\left(-s_0+s_1\left(A^2+\frac{1}{4}{B}^2\right)\right)=0,& (23)\end{array}$

which can equivalently be written by separating the real and imaginary parts as

$\begin{array}{ll}R-s_0+s_1\left(A^2+\frac{1}{4}{B}^2\right)=0,& (24)\end{array}$

$\begin{array}{ll}-{\omega }^2+\frac{1}{LC}=0.& (25)\end{array}$

Equations (22), (24), and (25) are then solved for A = Â = Q₀, $\omega =\widehat{\omega }=\frac{1}{\sqrt{LC}}$ and $B=\widehat{B}$ with $\widehat{B}$ such that

$\begin{array}{ll}{\widehat{B}}^2=4\left(\frac{s_0-R}{s_1}-Q_0^2\right)=4\left({\stackrel{-}{Q}}_0^2-Q_0^2\right).& (26)\end{array}$

Hence, for each index Q₀ such that $\left|Q_0\right|<{\stackrel{-}{Q}}_0$ there exists a PLC of the following form

$\begin{array}{ll}{q}_M(t)=Â+\widehat{B}\text{ cos }\widehat{\omega }t=Q_0+2\sqrt{{\stackrel{-}{Q^2}}_0-Q_0^2}\text{ cos }\frac{1}{\sqrt{LC}}t.& (27)\end{array}$

This implies that each manifold M(Q₀) contains a unique PLC for $\left|Q_0\right|<{\stackrel{-}{Q}}_0$. Since

$\begin{array}{ll}{i}_L(t)={D}{q}_M(t)=-\widehat{B}\widehat{\omega }\text{ sin }\widehat{\omega }t=-\frac{2}{\sqrt{LC}}\sqrt{{\stackrel{-}{Q^2}}_0-Q_0^2}\text{ sin }\frac{1}{\sqrt{LC}}t,& (28)\end{array}$

the trajectory described by the PLC in the (q_M, i_L)-plane lies on the following ellipse

$\begin{array}{ll}{\left(q_M-Q_0\right)}^2+LC{i}_L^2=4\left({\stackrel{-}{Q^2}}_0-Q_0^2\right).& (29)\end{array}$

Note that the ellipse is centered at the equilibrium point (Q₀, 0) and its size is maximum at Q₀ = 0 and decreases to zero as the index Q₀ tends to ${\stackrel{-}{Q}}_0$, thus collapsing to the equilibrium point. Moreover, since for $\left|Q_0\right|<{\stackrel{-}{Q}}_0$ the equilibrium point is unstable and taking into account that the system is two dimensional, we can conclude that the PLC is stable. In this regard, we mention that PLC stability can be assessed also via approximate HB tools, such as the Loeb criterion (see e.g., Atherton, 1975; Tesi et al., 1996).

Summing up, we have shown that in the input-less memristor circuit there coexist infinitely many attractors: a stable equilibrium point for each value of the manifold index Q₀ such that $\left|Q_0\right|>{\stackrel{-}{Q}}_0$ and a stable PLC for each value of the index Q₀ such that $\left|Q_0\right|<{\stackrel{-}{Q}}_0$. Figure 6 illustrates this multistability scenario for the following normalized values of the circuit parameters

$\begin{array}{ll}R=0.4,C=0.1,L=1.5,s_0=0.7,s_1=0.3,& (30)\end{array}$

from which it follows that ${\stackrel{-}{Q}}_0=1$.

FIGURE 6

Figure 6. Input-less circuit attractors: stable (green ⋆) equilibrium points and stable PLCs (solid green) given by Equations (27) and (28) as a function of the manifold index Q₀. The red circles denote the unstable equilibrium points.

In particular, it turns out that at $Q_0=\pm {\stackrel{-}{Q}}_0=\pm 1$ a typical behavior of (supercritical) Hopf bifurcation is observed. However, since it is obtained by varying the index Q₀ and thus the circuit initial conditions (for fixed circuit parameters), it may be referred to as a bifurcation without parameters (Corinto and Forti, 2017; Di Marco et al., 2018; Innocenti et al., 2019b; Ascoli et al., 2020a). As already discussed, it is expected that there exists a true limit cycle close to a PLC. Indeed, a more refined numerical analysis shows that for each |Q₀| < 1 the corresponding invariant manifold M(Q₀) has a unique limit cycle which attracts all the (non-trivial) trajectories. Moreover, the limit cycle is well-approximated by the PLC, as shown in Figure 7 for Q₀ = 0.

FIGURE 7

Figure 7. True limit cycle (solid dark) and PLC (solid red) in the (q_M, i_L)-plane for Q₀ = 0.

A quantitative comparison between the PLCs and the true limit cycles is provided by Figure 8 where the maximum and the minimum of both the true periodic solution q_M(t) and the PLC in Equation (27) are depicted as a function of the index Q₀.

FIGURE 8

Figure 8. Stable (green ⋆) and unstable (red ◦) equilibrium points; maximum (32) and minimum (31) of the PLCs (solid red) as a function of Q₀; maximum and minimum values of the numerically computed limit cycles (dotted dark points).

It can be readily checked that the minimum and the maximum of the first-order approximation are given by

$\begin{array}{ll}{\displaystyle \underset{t}{min}}{q}_M(t)=Q_0-2\sqrt{{\stackrel{-}{Q^2}}_0-Q_0^2}& (31)\end{array}$

and

$\begin{array}{ll}{\displaystyle \underset{t}{max}}{q}_M(t)=Q_0+2\sqrt{{\stackrel{-}{Q^2}}_0-Q_0^2},& (32)\end{array}$

respectively. We finally observe that similar results can be derived also for values of the circuit parameters different from those in Equation (30).

4. Modeling Neuron Dynamics via the Controlled Circuit

Section 3 makes it clear that in the input-less case the memristor circuit displays infinitely many attractors. Specifically, each planar invariant manifold M(Q₀) contains either a unique stable equilibrium point or a stable limit cycle. Moreover, an approximation of the limit cycle belonging to M(Q₀), with Q₀ such that $\left|Q_0\right|<{\stackrel{-}{Q}}_0$, is explicitly obtained in terms of the PLC (27).

In this section it will be shown how the coexistence of infinitely many stable equilibrium points and limit cycles, together with the knowledge of their dependence on the manifold index Q₀, which in the case of limit cycles is given in an approximate way by Equation (27), can be suitably exploited to make the memristor charge q_M responding as in Figure 2B to a pulse shaped input source I_s. To proceed, we consider the circuit state equations (3) by replacing v_C with the following state variable

$\begin{array}{ll}Q\doteq q_M+C{v}_C.& (33)\end{array}$

It turns out that (3) reduces to the third-order system

$\begin{array}{ll}\{\begin{array}{c}{D}Q(t)=I_s(t)\\ {D}{q}_M(t)=i_L\\ {D}{i}_L(t)=\frac{1}{LC}(Q-q_M)+\frac{s_0-R}{L}{i}_L-\frac{s_1}{L}{q}_M^2{i}_L\end{array}& (34)\end{array}$

where Q, q_M, and i_L are the new state variables. Suppose that at time t = t_i ≥ t₀ the current source I_s(t) displays a pulse of area Λ and width Δ > 0, i.e., I_s(t) is equal to zero for t ∈ [t₀, t_i] and t ∈ [t_i + Δ, +∞) and such that

$\begin{array}{ll}{\displaystyle {\int }_{t_i}^{t_i+\Delta }}{I}_s(t)dt=\Lambda .& (35)\end{array}$

From the first equation of Equation (34) it follows that

$\begin{array}{ll}Q(t)=Q_0^{(i)}+{\displaystyle {\int }_{t_i}^t}{I}_s(\sigma )d\sigma ,& (36)\end{array}$

which implies $Q(t)=Q_0^{(i)}$ for t ∈ [t₀, t_i] and $Q(t)=Q_0^{(i)}+\Lambda$ for t ∈ [t_i + Δ, +∞). For instance, in the case of a rectangular pulse we have $Q(t)=Q_0^{(i)}+(A/\Delta )(t-t_i)$.

By comparing the second and third equations of Equation (34) with those of Equation (7), it follows that the circuit dynamics lies onto $\mathit{\text{M}}(Q_0^{(i)})$ for t ∈ [t₀, t_i], it moves toward $\mathit{\text{M}}(Q_0^{(i)}+\Lambda )$ during the interval [t_i, t_i + Δ], it reaches $\mathit{\text{M}}(Q_0^{(i)}+\Lambda )$ at t = t_i + Δ, where it remains for all t ≥ t_i + Δ. Hence, the circuit has the property that pulse shaped current sources I_s make the dynamics moving from an initial manifold to any other one within a given time interval, by suitably choosing the pulse area and width.

We are interested in the case where the manifold $\mathit{\text{M}}(Q_0^{(i)})$ has a stable equilibrium point and in particular the index $Q_0^{(i)}$ is negative and hence $Q_0^{(i)}<-{\stackrel{-}{Q}}_0$. Moreover, we assume that at t = t_i the circuit state is at the equilibrium point, which means that $(Q(t_i),q_M(t_i),i_L(t_i))=(Q_0^{(i)},Q_0^{(i)},0)$. Let us now investigate the dynamics induced by a rectangular pulse of positive area Λ and width Δ on the circuit dynamics. It turns out that the final manifold $\mathit{\text{M}}(Q_0^{(i)}+\Lambda )$ still contains a stable equilibrium point if $\Lambda <\left|Q_0^{(i)}\right|-{\stackrel{-}{Q}}_0$, while it contains a stable limit cycle if $\left|Q_0^{(i)}\right|-{\stackrel{-}{Q}}_0<\Lambda <\left|Q_0^{(i)}\right|+{\stackrel{-}{Q}}_0$ and again a stable equilibrium point if $\Lambda >\left|Q_0^{(i)}\right|+{\stackrel{-}{Q}}_0$. Figure 9 illustrates the first two possible behaviors in the case of a rectangular pulse of height A/Δ with the circuit parameters as in (30).

FIGURE 9

Figure 9. Dynamics generated by applying at t_i = 10 a rectangular pulse of area Λ and width Δ = 1 via the current source I_s; the circuit parameters are as in Equation (30) and the initial conditions are q_M(0) = −2 and I_L(0) = 0, while Q(0) = −2. (A) For Λ = 0.5, the index Q(t) and the state variables q_M(t), i_L(t) converge to Q₀ = −1.5 and to the equilibrium point (−1.5, 0), respectively. (B) For Λ = 2.5, Q(t) converges to Q₀ = 0.5 and the state variables q_M(t), i_L(t) converge to the limit cycle belonging to the manifold M(0.5).

Specifically, after reaching $\mathit{\text{M}}(Q_0^{(i)}+\Lambda )$ at time t = t_i + Δ at the point (Q(t_i + Δ), q_M(t_i + Δ), i_L(t_i + Δ)) with $Q(t_i+\Delta )=Q_0^{(i)}+\Lambda$, the dynamics converges toward the attractor contained in the manifold, an equilibrium point in Figure 9A and a limit cycle in Figure 9B.

Clearly, the length of the transient behavior toward the attractor depends on the values of q_M(t_i + Δ) and i_L(t_i + Δ). The farther away is the point (q_M(t_i + Δ), i_L(t_i + Δ)) from the attractor lying onto $\mathit{\text{M}}(Q_0^{(i)}+\Lambda )$, the longer is the transient. Clearly, if (q_M(t_i + Δ), i_L(t_i + Δ)) belongs to the attractor we have no transient behavior. The values of q_M(t_i + Δ) and i_L(t_i + Δ) cannot be computed explicitly since this would require to integrate the second and third equations of (34) in the interval [t_i, t_i + Δ] with the initial condition $(q_M(t_i),i_L(t_i))=(Q_0^{(i)},0)$ and $Q(t)=Q_0^{(i)}+\Lambda (t-t_i)/\Delta$. However, by employing a Taylor series expansion for q_M(t) and i_L(t) for t ∈ [t_i, t_i + Δ], it turns out that

$\begin{array}{l}{q}_M(t_i+\Delta )=q_M(t_i)+{\alpha }_1\Delta +{\alpha }_2{\Delta }^2+\dots =Q_0^{(i)}+O(\Delta )\\ i_L(t_i+\Delta )=i_L(t_i)+{\beta }_1\Delta +{\beta }_2{\Delta }^2\dots ={\beta }_1\Delta +O({\Delta }^2),& (37)\end{array}$

where α_i and β_i, i = 1, 2 are constants. It is interesting to note that if the width Δ of the pulse is small, then for t ∈ [t_i, t_i + Δ] the charge q_M remains almost constant, while i_L varies from zero to a quantity proportional to Δ.

Summing up, for a suitable choice of the pulse area A and for sufficiently small width Δ, the pulse current I_s is able to steer, within the interval [t_i, t_i + Δ], the dynamics from the stable equilibrium point of $\mathit{\text{M}}(Q_0^{(i)})$ to the manifold $\mathit{\text{M}}(Q_0^{(i)}+A)$ containing a stable limit cycle. Moreover, during the time interval the charge q_M is almost constant and the current i_L remains close to zero.

We are now ready to show how it is possible to make the charge q_M display a behavior similar to that of Figure 2B in response to a pulse shaped input source I_s. Here, we are more interested in the dynamic response of the neuron than in its set and reset mechanisms. For these mechanisms we simply adopt the pulse timing of Figure 10 which is assumed to be generated via suitable logic gates.

FIGURE 10

Figure 10. Top: predicted behavior $q_M^0(t)$ of the memristor charge. Bottom: timing of the current pulses.

At time t = t_i the hardware detects a stimulus and generates a current pulse I_s with area A and width Δ (set mechanism), which moves the dynamics away from the stable equilibrium point (resting point) in $\mathit{\text{M}}(Q_0^{(i)})$; at time t = t_i + Δ the dynamics reaches $\mathit{\text{M}}(Q_0^{(i)}+\Lambda )$ and the memristor charge q_M starts displaying a behavior similar to that reported in Figure 2B; at time t = t_i + Δ + T an inverse current pulse is generated in order to make the dynamics come back to the stable equilibrium point in $\mathit{\text{M}}(Q_0^{(i)})$ (reset mechanism).

Hence, the problem is how the parameters $Q_0^{(i)}$, Λ, Δ, and T should be designed in order to obtain the sought behavior of q_M. Let Δ be chosen enough small to ensure, as discussed above, that in [t_i, t_i + Δ] the charge q_M is almost constant and the current i_L remains close to zero. Choosing a small Δ implies that when at t = t_i + Δ the dynamics reaches the manifold $\mathit{\text{M}}(Q_0^{(i)}+\Lambda )$ we have that q_M and i_L are quite close to $Q_0^{(i)}$ and zero, respectively. This manifold contains a stable limit cycle which can be in general computed only numerically. However, in section 3 it has been shown that the limit cycle can be approximated by a PLC. Clearly, the current expression of the PLC on $\mathit{\text{M}}(Q_0^{(i)}+\Lambda )$ is obtained by putting $Q_0=Q_0^{(i)}+\Lambda$ in Equation (27).

Now, observe that the minimum and the maximum values of q_M, expressed in Equations (31) and (32), respectively, are achieved once i_L is equal to zero. Hence, since for sufficiently small Δ we have that $(q_M(t_i+\Delta ),i_L(t_i+\Delta ))\approx (Q^{(i)},0)$, it is possible to set (q_M(t_i + Δ), i_L(t_i + Δ)) as the point of the manifold $\mathit{\text{M}}(Q_0^{(i)}+\Lambda )$ where the PLC has its minimum value (31). Consequently, for t ∈ [t_i + Δ, t_i + Δ + T] with T being the PLC period, the PLC provides the typical harmonic oscillator over one period, i.e., it evolves from the minimum to the maximum coming back again to the minimum. At the end of the period, the reset mechanism generates a current pulse of area −A and width Δ, which makes the dynamics come back to the stable equilibrium point of $\mathit{\text{M}}(Q_0^{(i)})$. Hence, the resulting behavior of the memristor charge $q_M^0(t)$ can be expressed as

$\begin{array}{ll}{q}_M^0(t)=\{\begin{array}{ll}{Q}_0^{(i)}\hfill & \hfill t\in [t_0,t_i+\Delta ]\hfill \\ Q_0^{(i)}+\Lambda +2\sqrt{{\overline{Q^2}}_0-{(Q_0^{(i)}+\Lambda )}^2}\cos \frac{1}{\sqrt{LC}}(t-t_i-\Delta )+\pi )\hfill & \hfill t\in (t_i+\Delta ,t_i+\Delta +T)\hfill \\ Q_0^{(i)}\hfill & \hfill t\in (t_i+\Delta +T,+\infty )\hfill \end{array}& (38)\end{array}$

and is depicted in Figure 10. Note that the parameter T has been chosen as the period of the PLC, i.e.:

$\begin{array}{ll}T=2\pi \sqrt{LC}.& (39)\end{array}$

Hereafter, we denote by $i_L^0(t)$ the time-derivative of $q_M^0(t)$, i.e.:

$\begin{array}{ll}{i}_L^0(t)={D}{q}_M^0(t)\forall t\in [t_0,+\infty ).& (40)\end{array}$

Now, to set (q_M(t_i + Δ), i_L(t_i+Δ)) as the point of the manifold $\mathit{\text{M}}(Q_0^{(i)}+\Lambda )$ where the PLC has its minimum value (31), the following condition must be satisfied

$\begin{array}{ll}{Q}_0^{(i)}=Q_0-2\sqrt{{\stackrel{-}{Q}}_0^2-Q_0^2}& (41)\end{array}$

with $Q_0=Q_0^{(i)}+\Lambda$. It is not difficult to verify that there are many solutions $Q_{(i)}^0$ and Q₀ of (41) such that $Q_0^{(i)}<-{\stackrel{-}{Q}}_0$ and $\left|Q_0\right|<{\stackrel{-}{Q}}_0$. We are interested in the one with the minimum value of $Q_0^{(i)}$ because this choice ensures that the stable equilibrium point on the initial manifold $\mathit{\text{M}}(Q_0^{(i)})$ is at the maximum distance from the invariant manifold $\mathit{\text{M}}({\stackrel{-}{Q}}_0)$ where the Hopf bifurcation happens, thus guaranteeing a larger robustness margin against spurious pulses of small area. The minimum value of $Q_0^{(i)}$ can be computed in closed form by minimizing the right hand side of Equation (41). It turns out that

$\begin{array}{ll}{Q}_0\doteq Q_0^{*}=-\frac{1}{\sqrt{5}}{\stackrel{-}{Q}}_0,& (42)\end{array}$

and since

$\begin{array}{l}{Q}_0^{*}-2\sqrt{{\stackrel{-}{Q}}_0^2-{(Q_0^{*})}^2}=-\sqrt{5}{\stackrel{-}{Q}}_0,\end{array}$

from Equations (41) and (9) we have

$\begin{array}{ll}{Q}_0^{(i)}=-\sqrt{5}\sqrt{\frac{s_0-R}{s_1}}.& (43)\end{array}$

Moreover, since Q₀ in Equation (41) stands for $Q_0^{(i)}+A$, from Equation (42) we finally get

$\begin{array}{ll}\Lambda =\frac{4}{\sqrt{5}}\sqrt{\frac{s_0-R}{s_1}}.& (44)\end{array}$

Finally, the pulse width Δ should be practically chosen to be some percent of the period T. Hence, we can select it according to the relation

$\begin{array}{ll}\Delta =2\sigma \pi \sqrt{LC},& (45)\end{array}$

where σ lies in the interval [0.01, 0.05].

Summing up, to obtain $q_M^0(t)$ in Equation (38), the parameters T, $Q_0^{(i)}$, Λ, and Δ can be designed according to the corresponding expressions given in Equations (39), (43), (44), and (45), respectively. Notably, these relations depend explicitly on the circuit parameters R, L, C, s₀, and s₁. Hence, by suitably adapting these parameters it is possible to vary the features of $q_M^0(t)$.

On the other hand, the formulas in Equations (39), (43), (44), and (45) have been derived under two approximations. The first one concerns the assumption that at time t = t_i + Δ the values of q_M and i_L are quite close to $Q_0^{(i)}$ and zero, respectively. However, according to Equation (37) by lowering the value of the pulse width this assumption can be suitably satisfied. The second approximation is that the formulas have been obtained by using the PLCs instead of the true limit cycles. Indeed, as shown in section 3.2, the true limit cycles contain higher order harmonics terms which can generate some distortion in both the period and the minimum and maximum amplitudes over the period, with respect to those analytically computed for the PLCs. However, these higher order harmonics can be filtered out by suitably adjusting the circuit parameters. Moreover, the fact that the input-less memristor circuit contains a bundle of limit cycles makes the formulas, and hence the parameter design procedure, sufficiently robust with respect to the use of PLC in their derivation. Said another way, by slightly varying the parameters T, $Q_0^{(i)}$, A, and Δ from the nominal values provided by Equations (39), (43), (44), and (45), respectively, the distortion effect could be reduced.

To illustrate the design procedure let us consider the circuit parameters in Equation (30). From Equations (39), (43), and (44), we have that

$\begin{array}{ll}T=2.4334,Q_0^{(i)}=-2.2361,\Lambda =1.7889,& (46)\end{array}$

and choosing σ = 0.02 in Equation (45) we obtain

$\begin{array}{ll}\Delta =0.0487.& (47)\end{array}$

The corresponding predicted behavior $q_M^0(t)$ in Equation (38) is reported Figure 11A together with the true behavior of q_M(t) which is obtained by solving numerically (34). The state space trajectories corresponding to the true solution (Q(t), q_M(t), i_L(t)) of the circuit and the predicted solution $(Q(t),q_M^0(t),i_L^0(t))$) are depicted in Figure 11B.

FIGURE 11

Figure 11. The set pulse is applied at t_i = 2 with area A = 1.7889 and width Δ = 0.0487, the reset pulse is applied at t_i + Δ + T = 4.4821 with area Λ = −1.7889 and width Δ = 0.0487. The initial conditions are Q(0) = −2.2361, q_M(0) = −2.2361 and i_L(0) = 0. (A) Time behaviors of q_M (blu) and $q_M^0$ (red). (B) Trajectories in the state space generated by (Q(t), q_M(t), i_L(t)) (blu) and $(Q(t),q_M^0(t),i_L^0(t))$ (red).

Note that the charge q_M(t) has a dynamical behavior that is quite similar to the one depicted in Figure 2B, thus clearly showing that the memristor circuit can be used for mimicking some features of neuron dynamics.

This dynamical similarity is confirmed also in other scenarios. For instance, suppose that in the pulse timing the reset pulse is activated at t = t_i + Δ + nT, with n being a positive integer. Clearly, in this case $q_M^0(t)$ completes n periods within the interval t_i + Δ, t_i + Δ + nT, and, therefore, q_M(t) is expected to generate a burst of n spikes. Figure 12 provides the numerically computed behavior of q_M(t) for n = 8 over a time frame covering three bursts.

FIGURE 12

Figure 12. The initial conditions are as in Figure 11, the set pulses are activated at t_i ∈ {10, 60, 110} and the reset pulses are applied 19.5164 time units later. (A) Time behavior of q_M for n = 8. (B) Time behavior of I_s.

5. Conclusions

In this paper, a memristor circuit composed of a resistor, an inductor, a capacitor, an ideal charge-controlled memristor and an independent current source as input is considered. It is first shown that in the input-less case the circuit enjoys the foliation property of the state space, i.e., it contains infinitely many planar invariant manifolds which are parameterized by a scalar index depending on the circuit initial conditions. Each manifold contains an attractor which can be either a stable equilibrium point or a stable limit cycle, depending on the value of the manifold index. Moreover, a first-order periodic approximation is obtained in an analytic way for each limit cycle via the Describing Function (DF) technique, a classical tool within the Harmonic Balance (HB) context.

Then, it is shown that the memristor charge can mimic a simplified model of a neuron response when an external independent pulse-programmed current source is introduced in the circuit. Specifically, the sought dynamics of the memristor charge is generated via the concatenation of convergent and oscillatory behaviors, which are obtained by switching between stable equilibrium points and limit cycles via a suitable design of the pulse timing of the current source. Some relationships between the pulse and the circuit parameters are also devised exploiting the knowledge of the first-order periodic approximation of the limit cycles.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available, without undue reservation, on request to the corresponding author with appropriate motivation.

Author Contributions

All authors contributed to the conception and design of the study, wrote the manuscript, read, and approved the submitted version.

Funding

This work was supported by the project Analogue COmputing with Dynamic Switching Memristor Oscillators: Theory, Devices and Applications (COSMO), Grant PRIN 2017LSCR4K 002, funded by the Italian Ministry of Education, University and Research (MIUR).

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. ^Throughout the paper, ${D}f(\stackrel{-}{t})$ denotes the derivative of the function f(t) at $t=\stackrel{-}{t}$.

2. ^Consider a linear time-invariant finite-dimensional causal system described by the input-output relation $y(t)=L({D})u(t)$ where the operator $L({D})$ has the following general form

$L({D})=\frac{b_{n-1}{{D}}^{n-1}+\dots +b_1{D}+b_0}{{{D}}^n+a_{n-1}{{D}}^{n-1}+\dots +a_1{D}+a_0}$

with a_i, b_i, i = 0, …, n − 1 being constants. Let u(t) be an exponential signal, i.e., u(t) = Ce^σt, C ∈ ℝ, σ ∈ ℂ, and assume that L(σ) is finite, i.e., L(σ) ≠ ∞. Then, for suitable initial conditions of the system, we have y(t) = CL(σ)e^σt, i.e., the output of the system has the same exponential form of the input.

References

Adhikari, S. P., Yang, C., Kim, H., and Chua, L. O. (2012). Memristor bridge synapse-based neural network and its learning. IEEE Trans. Neural Netw. Learn. Syst. 23, 1426–1435. doi: 10.1109/TNNLS.2012.2204770

PubMed Abstract | CrossRef Full Text | Google Scholar

Ascoli, A., Messaris, I., Tetzlaff, R., and Chua, L. O. (2020a). Theoretical foundations of memristor cellular nonlinear networks: stability analysis with dynamic memristors. IEEE Trans. Circ. Syst. I Reg. Pap. 67, 1389–1401. doi: 10.1109/TCSI.2019.2957813

CrossRef Full Text | Google Scholar

Ascoli, A., Tetzlaff, R., Kang, S. M., and Chua, L. O. (2020b). Theoretical foundations of memristor cellular nonlinear networks: A DRM₂-based method to design memcomputers with dynamic memristors. IEEE Trans. Circ. Syst. I Reg. Pap. 67, 2753–2766. doi: 10.1109/TCSI.2020.2978460

CrossRef Full Text | Google Scholar

Atherton, D. P. (1975). Nonlinear Control Engineering. London: Van Nostrand Reinhold.

Google Scholar

Babacan, Y., Kaçar, F., and Gürkan, K. (2016). A spiking and bursting neuron circuit based on memristor. Neurocomputing 203, 86–91. doi: 10.1016/j.neucom.2016.03.060

CrossRef Full Text | Google Scholar

Basso, M., Genesio, R., and Tesi, A. (1997). A frequency method for predicting limit cycle bifurcations. Nonlinear Dyn. 13, 339–360. doi: 10.1023/A:1008298205786

CrossRef Full Text | Google Scholar

Bonani, F., and Gilli, M. (1999). Analysis of stability and bifurcations of limit cycles in Chua's circuit through the harmonic-balance approach. IEEE Trans. Circ. Syst. I Fundament. Theory Appl. 46, 881–890. doi: 10.1109/81.780370

CrossRef Full Text | Google Scholar

Bonomi, F., Milito, R., Zhu, J., and Addepalli, S. (2012). “Fog computing and its role in the internet of things,” in Proceedings of the First Edition of the MCC Workshop on Mobile Cloud Computing (Helsinki: ACM), 13–16. doi: 10.1145/2342509.2342513

CrossRef Full Text | Google Scholar

Chen, M., Sun, M., Bao, B., Wu, H., Xu, Q., and Wang, J. (2018). Controlling extreme multistability of memristor emulator-based dynamical circuit in flux-charge domain. Nonlinear Dyn. 91, 1395–1412. doi: 10.1007/s11071-017-3952-9

CrossRef Full Text | Google Scholar

Chen, Z., Zheng, M., Friesen, W. O., and Iwasaki, T. (2008). Multivariable harmonic balance analysis of the neuronal oscillator for leech swimming. J. Comput. Neurosci. 25:583. doi: 10.1007/s10827-008-0105-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Chua, L., Sbitnev, V., and Kim, H. (2012). Hodgkin-Huxley axon is made of memristors. Int. J. Bifurc. Chaos 22:1230011. doi: 10.1142/S021812741230011X

CrossRef Full Text | Google Scholar

Chua, L. O. (1971). Memristor-the missing circuit element. IEEE Trans. Circ. Theory 18, 507–519. doi: 10.1109/TCT.1971.1083337

CrossRef Full Text | Google Scholar

Corinto, F., Ascoli, A., and Gilli, M. (2011). Nonlinear dynamics of memristor oscillators. IEEE Trans. Circ. Syst. I Reg. Pap. 58, 1323–1336. doi: 10.1109/TCSI.2010.2097731

CrossRef Full Text | Google Scholar

Corinto, F., Di Marco, M., Forti, M., and Chua, L. (2019). Nonlinear networks with mem-elements: complex dynamics via flux-charge analysis method. IEEE Trans. Cybern. 1–14. doi: 10.1109/TCYB.2019.2904903

PubMed Abstract | CrossRef Full Text | Google Scholar

Corinto, F., and Forti, M. (2016). Memristor circuits: flux-charge analysis method. IEEE Trans. Circ. Syst. I Reg. Pap. 63, 1997–2009. doi: 10.1109/TCSI.2016.2590948

CrossRef Full Text | Google Scholar

Corinto, F., and Forti, M. (2017). Memristor circuits: bifurcations without parameters. IEEE Trans. Circ. Syst. I Reg. Pap. 64, 1540–1551. doi: 10.1109/TCSI.2016.2642112

CrossRef Full Text | Google Scholar

Corinto, F., and Forti, M. (2018). Memristor circuits: pulse programming via invariant manifolds. IEEE Trans. Circ. Syst. I Reg. Pap. 65, 1327–1339. doi: 10.1109/TCSI.2017.2740999

CrossRef Full Text | Google Scholar

Di Marco, M., Forti, M., Innocenti, G., and Tesi, A. (2018). Harmonic balance method to analyze bifurcations in memristor oscillatory circuits. Int. J. Circ. Theory Appl. 46, 66–83. doi: 10.1002/cta.2414

CrossRef Full Text | Google Scholar

Di Marco, M., Forti, M., Innocenti, G., and Tesi, A. (2020a). “Input design for controlling dynamics in a second-order memristive circuit,” in 2020 European Conference on Circuit Theory and Design (ECCTD) (Sofia), 1–4. doi: 10.1109/ECCTD49232.2020.9218368

CrossRef Full Text | Google Scholar

Di Marco, M., Forti, M., Innocenti, G., and Tesi, A. (2021). “Transient control in targeting multistable dynamics of a memristor circuit,” in 2021 IEEE International Symposium on Circuits and Systems (ISCAS 2021) (Daegu: IEEE). doi: 10.1109/ISCAS51556.2021.9401351

CrossRef Full Text | Google Scholar

Di Marco, M., Forti, M., Innocenti, G., Tesi, A., and Corinto, F. (2020b). “Targeting multistable dynamics in a second-order memristor circuit,” in 2020 IEEE International Symposium on Circuits and Systems (ISCAS) (Seville), 1–5. doi: 10.1109/ISCAS45731.2020.9181111

CrossRef Full Text | Google Scholar

Di Marco, M., Forti, M., and Pancioni, L. (2016). Complete stability of feedback CNNs with dynamic memristors and second-order cells. Int. J. Circ. Theory Appl. 44, 1959–1981. doi: 10.1002/cta.2205

CrossRef Full Text | Google Scholar

Di Marco, M., Forti, M., and Pancioni, L. (2017). Convergence and multistability of nonsymmetric cellular neural networks with memristors. IEEE Trans. Cybern. 47, 2970–2983. doi: 10.1109/TCYB.2016.2586115

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Di Marco, M., Forti, M., and Tesi, A. (2003). Harmonic balance approach to predict period-doubling bifurcations in nearly-symmetric neural networks. J. Circ. Syst. Comput. 12, 435–460. doi: 10.1142/S0218126603000969

CrossRef Full Text | Google Scholar

Di Marco, M., Innocenti, G., Forti, M., and Tesi, A. (2019). “Control design for targeting dynamics of memristor Murali-Lakshmanan-Chua circuit,” in IEEE ECC 2019 (Naples), 4332–4337. doi: 10.23919/ECC.2019.8795813

CrossRef Full Text | Google Scholar

Du, C., Cai, F., Zidan, M. A., Ma, W., Lee, S. H., and Lu, W. D. (2017). Reservoir computing using dynamic memristors for temporal information processing. Nat. Comm. 8:2204. doi: 10.1038/s41467-017-02337-y

PubMed Abstract | CrossRef Full Text

Gelb, A., and Vander Velde, W. E. (1968). Multiple-Input Describing Functions and Nonlinear System Design. New York, NY: McGraw-Hill Book Co.

Google Scholar

Genesio, R., and Tesi, A. (1992). Harmonic balance methods for the analysis of chaotic dynamics in nonlinear systems. Automatica 28, 531–548. doi: 10.1016/0005-1098(92)90177-H

CrossRef Full Text | Google Scholar

Hong, Q., Zhao, L., and Wang, X. (2019). Novel circuit designs of memristor synapse and neuron. Neurocomputing 330, 11–16. doi: 10.1016/j.neucom.2018.11.043

CrossRef Full Text | Google Scholar

Hu, M., Chen, Y., Yang, J. J., Wang, Y., and Li, H. H. (2016). A compact memristor-based dynamic synapse for spiking neural networks. IEEE Trans. Comput. Aided Design Integr. Circ. Syst. 36, 1353–1366. doi: 10.1109/TCAD.2016.2618866

CrossRef Full Text | Google Scholar

Ielmini, D., and Pedretti, G. (2020). Device and circuit architectures for in-memory computing. Adv. Intell. Syst. 2:2000040. doi: 10.1002/aisy.202000040

CrossRef Full Text | Google Scholar

Ielmini, D., and Wong, H.-S. P. (2018). In-memory computing with resistive switching devices. Nat. Electron. 1:333. doi: 10.1038/s41928-018-0092-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Innocenti, G., Di Marco, M., Forti, M., and Tesi, A. (2019a). “A controlled Murali-Lakshmanan-Chua memristor circuit to mimic neuron dynamics,” in 2019 IEEE 58th Conference on Decision and Control (CDC) (Nice: IEEE), 3423–3428. doi: 10.1109/CDC40024.2019.9029330

CrossRef Full Text | Google Scholar

Innocenti, G., Di Marco, M., Forti, M., and Tesi, A. (2019b). Prediction of period doubling bifurcations in harmonically forced memristor circuits. Nonlinear Dyn. 96, 1169–1190. doi: 10.1007/s11071-019-04847-4

CrossRef Full Text | Google Scholar

Innocenti, G., Di Marco, M., Tesi, A., and Forti, M. (2020). Input-output characterization of the dynamical properties of circuits with a memelement. Int. J. Bifurc. Chaos 30:2050110. doi: 10.1142/S0218127420501102

CrossRef Full Text | Google Scholar

Innocenti, G., Tesi, A., and Genesio, R. (2010). Complex behaviour analysis in quadratic jerk systems via frequency domain Hopf bifurcation. Int. J. Bifurc. Chaos 20, 657–667. doi: 10.1142/S0218127410025946

CrossRef Full Text | Google Scholar

Izhikevich, E. M. (2000). Neural excitability, spiking and bursting. Int. J. Bifurc. Chaos 10, 1171–1266. doi: 10.1142/S0218127400000840

CrossRef Full Text | Google Scholar

Izhikevich, E. M. (2007). Dynamical Systems in Neuroscience. Cambridge, MA: MIT Press. doi: 10.7551/mitpress/2526.001.0001

CrossRef Full Text | Google Scholar

Jo, S. H., Chang, T., Ebong, I., Bhadviya, B. B., Mazumder, P., and Lu, W. (2010). Nanoscale memristor device as synapse in neuromorphic systems. Nano Lett. 10, 1297–1301. doi: 10.1021/nl904092h

PubMed Abstract | CrossRef Full Text | Google Scholar

Khalil, H. K. (2002). Nonlinear Systems. Upper Saddle River, NJ: Prentice Hall.

Kim, S., Du, C., Sheridan, P., Ma, W., Choi, S., and Lu, W. D. (2015). Experimental demonstration of a second-order memristor and its ability to biorealistically implement synaptic plasticity. Nano Lett. 15, 2203–2211. doi: 10.1021/acs.nanolett.5b00697

PubMed Abstract | CrossRef Full Text | Google Scholar

Krestinskaya, O., James, A. P., and Chua, L. O. (2020). Neuromemristive circuits for edge computing: a review. IEEE Trans. Neural Netw. Learn. Syst. 31, 4–23. doi: 10.1109/TNNLS.2019.2899262

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, S., Strachan, J. P., and Williams, R. S. (2017). Chaotic dynamics in nanoscale NbO₂ Mott memristors for analogue computing. Nature 548:318. doi: 10.1038/nature23307

CrossRef Full Text | Google Scholar

Matsuoka, K. (2011). Analysis of a neural oscillator. Biol. Cybernet. 104, 297–304. doi: 10.1007/s00422-011-0432-z

CrossRef Full Text | Google Scholar

Mees, A. I. (1981). Dynamics of Feedback Systems. New York, NY: Wiley.

Google Scholar

Nakada, K. (2019). An ReRAM-based neuron device for neuromorphic pulse coding. Jpn. J. Appl. Phys. 58:030904. doi: 10.7567/1347-4065/aafb60

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Piccardi, C. (1994). Bifurcations of limit cycles in periodically forced nonlinear systems: the harmonic balance approach. IEEE Trans. Circ. Syst. I 41, 315–320. doi: 10.1109/81.285687

CrossRef Full Text | Google Scholar

Pisarchik, A. N., and Feudel, U. (2014). Control of multistability. Phys. Rep. 540, 167–218. doi: 10.1016/j.physrep.2014.02.007

CrossRef Full Text | Google Scholar

Satyanarayanan, M. (2017). The emergence of edge computing. Computer 50, 30–39. doi: 10.1109/MC.2017.9

CrossRef Full Text | Google Scholar

Tesi, A., Abed, E. H., Genesio, R., and Wang, H. O. (1996). Harmonic balance analysis of period-doubling bifurcations with implications for control of nonlinear dynamics. Automatica 32, 1255–1271. doi: 10.1016/0005-1098(96)00065-9

CrossRef Full Text | Google Scholar

Thomas, A. (2013). Memristor-based neural networks. J. Phys. D Appl. Phys. 46:093001. doi: 10.1088/0022-3727/46/9/093001

CrossRef Full Text | Google Scholar

Usha, K., and Subha, P. (2019). Hindmarsh-Rose neuron model with memristors. Biosystems 178, 1–9. doi: 10.1016/j.biosystems.2019.01.005

CrossRef Full Text | Google Scholar

Waldrop, M. M. (2016). The chips are down for Moore's law. Nat. News 530:144. doi: 10.1038/530144a

CrossRef Full Text | Google Scholar

Wang, N., Zhang, G., and Bao, H. (2019). Bursting oscillations and coexisting attractors in a simple memristor-capacitor-based chaotic circuit. Nonlinear Dyn. 97, 1477–1494. doi: 10.1007/s11071-019-05067-6

CrossRef Full Text | Google Scholar

Wang, Z., Ambrogio, S., Balatti, S., and Ielmini, D. (2015). A 2-transistor/1-resistor artificial synapse capable of communication and stochastic learning in neuromorphic systems. Front. Neurosci. 8:438. doi: 10.3389/fnins.2014.00438

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Xu, Q., Lin, Y., Bao, B., and Chen, M. (2016). Multiple attractors in a non-ideal active voltage-controlled memristor based Chua's circuit. Chaos Solitons Fractals 83, 186–200. doi: 10.1016/j.chaos.2015.12.007

CrossRef Full Text | Google Scholar

Yang, J. J., Strukov, D. B., and Stewart, D. R. (2013). Memristive devices for computing. Nat. Nanotechnol. 8, 13–24. doi: 10.1038/nnano.2012.240

CrossRef Full Text | Google Scholar

Yuan, F., Wang, G., Jin, P., Wang, X., and Ma, G. (2016). Chaos in a meminductor-based circuit. Int. J. Bifurc. Chaos 26:1650130. doi: 10.1142/S0218127416501303

CrossRef Full Text | Google Scholar

Zidan, M. A., Strachan, J. P., and Lu, W. D. (2018). The future of electronics based on memristive systems. Nat. Electron. 1:22. doi: 10.1038/s41928-017-0006-8

CrossRef Full Text | Google Scholar