Embedding the Dynamics of Forced Nonlinear Systems in Multistable Memristor Circuits

A well-known feature of memristors is that they makes the circuit dynamics much richer than that generated by classical $RLC$ circuits containing nonlinear resistors. In the case of circuits with ideal memristors, such a multistability property, i.e., the coexistence of many different attractors for a fixed set of parameters, is connected to the fact that the state space is composed of a continuum of invariant manifolds where either convergent or oscillatory and more complex behaviors can be displayed. In this paper we investigate the possibility of designing memristor circuits where known attractors are embedded into the invariant manifolds. We consider a class of forced nonlinear systems containing several systems which are known to display complex dynamics, and we investigate under which conditions the dynamics of any given system of the class can be reproduced by a circuit composed of a two-terminal (one port) element connected to a flux-controlled memristor. It is shown that an input-less circuit is capable to replicate the system attractors generated by varying the constant forcing input, once the parameters of the two-terminal element and the memristor nonlinear characteristic are suitably selected. Indeed, there is a one-one correspondence between the dynamics generated by the nonlinear system for a constant value of the input and that displayed on one of the invariant manifolds of the input-less memristor circuit. Some extensions concerning the case of non-constant forcing terms and the use of charge-controlled memristors are also provided. The results are illustrated via FitzHugh-Nagumo model and Duffing oscillator.


I. INTRODUCTION
R ECENT years have witnessed a widespread interest in exploring possible ways to alleviate some emerging limitations of digital Von Neumann computing systems [1], [2]. Within this context, in-memory computing is seen as a promising approach to overcome the memory bottleneck related to the strongly increasing amount of data to be processed in solving computational tasks. In-memory computing is based on new nanoscale devices, such as memristors, possessing unconventional features which make it possible to foreseen the development of novel analogue and parallel (neuromorphic) computing schemes [3], [4].
The ideal memristor has been introduced in the seminal 1971 paper by Chua [5] as the fourth basic passive circuit element in addition to the resistor, capacitor and inductor. Its constitutive relation is a nonlinear function linking charge to flux (resp., flux to charge) for a flux-controlled (resp., charge-controlled) memristor. In the voltage-current domain an ideal memristor satisfies a state-dependent Ohm's law where the state variable is either the flux, for a flux-controlled memristor, or the charge, for a charge-controlled memristor. Later on, a wider class of memristive systems has been proposed by Chua and Kang [6]. These have been further subdivided in generic and extended memristors, depending upon the complexity of the constitutive relation involving voltage, current and additional memristor internal state variables. We refer the reader to [7] for the nomenclature and a discussion of the hierarchy and genealogy of memristor models used in the literature. Generic and extended memristors are of practical importance in that they can be used to better model some classes of real memristive devices in nanotechnology with respect to ideal memristors [8], [9], [10]. On the other hand, ideal memristors are of great interest in nonlinear circuit theory [11], [12], [13], [14]. Moreover, several articles in the literature discuss physical components or systems whose dynamic behavior can be approximated by that of ideal memristors [15], [16], [17]. In particular, [18] and [19] show that in practical ranges of temperature and thickness of the amorphous region, mushroom-type phase-change memory (PCM) and ReRAM devices obey a flux-charge model as that describing ideal memristors. It is also worth to note that numerous analog or digital techniques are available to implement ideal memristors. We refer the reader to [13] and [20] for a detailed account of the current state of the art. Such techniques are based on using current-mode building blocks as second-generation Current Conveyor (CCII), electronically tunable CCII, operational transconductance amplifier (OTA), current feedback operational amplifier (CFOA), differential difference current conveyor, current conveyor transconductance amplifier (CCTA), and differential voltage CCTA. Also, relatively simple structures based on MOS transistors have been proposed (see, e.g., [20] and references therein). Wavedigital emulators of memristors have been devised in [21].
From a dynamic point of view, the key feature of memristors, either ideal or extended and real physical ones, is that their presence makes the circuits capable to generate a large This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ variety of attractors for a fixed set of circuit parameters [12], [22], [23], [24], [25]. This important property is referred to as (extreme) multistability and it is not enjoyed by traditional nonlinear RLC circuits which generally display only a limited number of attractors for any fixed set of circuit parameters. For ideal memristors, it has been shown that multistability is connected to the fact that the state space of memristor circuits is decomposed into a continuum of invariant manifolds [11]. Since on each manifold either convergent or oscillatory and more complex behaviors can be displayed, the coexistence of infinitely many attractors is a natural scenario for circuits with ideal memristors [26], [27], [28], [29], [30], [31], [32]. Moreover, each invariant manifold can be uniquely identified via some parameters referred to as manifold indexes, whose values depend on the initial conditions of the circuit. By suitably setting these indexes, it is possible to select a desired manifold and also to switch the circuit dynamics between different invariant manifolds and the attractors therein contained [33], [34], [35], [36], [37]. In particular, [37] shows how a shaped voltage/current source can be pulse programmed in order to steer the dynamics from one invariant manifold to another.
The multistability feature of memristors has been thoroughly investigated, also in connection with some classic problems (see, e.g., [38] and [39]). However, several interesting issues still deserve to be analyzed. One concerns how many and which type of attractors can be embedded in a memristor circuit, an issue which is certainly appealing from a theoretical point of view since it leads to a better understanding of the richness of the memristor circuits dynamics. It can also have some practical implications, as for instance in the realization of reservoir computing systems [40]. The problem of embedding known dynamics in some system is indeed a classical problem. For instance, a fundamental result of Smale ensures that any dynamic system of order n − 1 can be embedded in a competitive system of order n (see, e.g., [41] at pp. 344-345). However, this competitive system has not the structure of a circuit or a neural network. Here, we are interested in analyzing how many attractors can be embedded in a system which is implementable via a memristor circuit. Another issue is connected with the role of the input to generate attractors in forced nonlinear systems. For instance, it is well-known that in the FitzHugh-Nagumo model [44], [45] either convergent or oscillatory behaviors are displayed depending on the constant value of the input, i.e., the injected current.
In this paper we investigate how the multistability feature of circuits with an ideal memristor can be exploited to reproduce all the dynamics displayed by a forced nonlinear system. If the forcing input is constant, as in the case of the FitzHugh-Nagumo model, then we consider an input-less memristor circuit whose manifold index is related to the value of the input. If the input is non-constant, as for instance in the case of the Duffing system, then some independent voltage/current sources are introduced in the circuit and programmed according to the index variations. Specifically, the considered class of forced nonlinear systems contains several systems which are known to display complex dynamics, and it is described in Section II together with the memristor circuit to be designed. Such a circuit is assumed to be composed of a finite-dimensional causal linear time-invariant two-terminal (one port) element connected to a flux-controlled ideal memristor. Section III investigates under which conditions an input-less memristor circuit is capable to exactly reproduce the dynamics generated when the nonlinear system is forced by a constant input. Specifically, we look for conditions on the circuit structure which ensure the existence of a one-toone correspondence between the dynamics of the nonlinear system for a given value of the constant forcing input and that displayed on one of the invariant manifolds. Moreover, a procedure to synthesize the two-terminal element is discussed and its application to the FitzHugh-Nagumo model is illustrated in Section IV. Section V considers the extension to the case when the forcing input of the nonlinear system is no longer constant, by showing that the system dynamics can be still replicated once the two-terminal element is equipped with suitable programmed voltage/current sources. Finally, Section VI considers the case when the flux-controlled memristor is replaced by a charge-controlled one, while some concluding remarks are given in Section VII.

II. PRELIMINARIES AND PROBLEM FORMULATION
In this paper we consider the class of nonlinear systems which are described by the following state space representation : where x ∈ R n is the state vector, u ∈ R is a scalar input, y ∈ R is a scalar output, A ∈ R n×n , B ∈ R n×1 , C ∈ R 1×n , D ∈ R n×1 . We assume that the nonlinear function f : R → R is such that f (0) = 0 and it is sufficiently smooth to ensure existence and uniqueness of the solutions of . The following assumptions are enforced on the matrices A, B, and C. Assumption 1: The pair (A, B) is controllable and the pair (A, C) is observable.
Also, we initially assume the matrix A is non-singular and the input u is constant over time, i.e., We observe that the structure of is enjoyed by many systems displaying either convergent or oscillatory and more complex behaviors, such as the FitzHugh-Nagumo model, the Van der Pol and Duffing oscillators, the Chua and Murali-Lakshmanan-Chua circuits. In the sequel we denote by x 0 (t) the solution of (1) for t ≥ t 0 with initial conditions x 0 (t 0 ) = x 0 and constant input U 0 . Accordingly, we refer to S as the set of all the solutions of generated by varying x 0 ∈ R n and U 0 ∈ R. We find it convenient to write this set as where S U 0 is the set of all the solutions x 0 (t) obtained for the fixed input u(t) = U 0 by varying the initial conditions Our aim is to show that the set S can be reproduced via the dynamics of the circuit depicted in Fig. 1, which is composed of a two-terminal (one port) element L, with voltage v L and current i L , and a memristor MR, with voltage v M and current i M . The two-terminal element L is either a passive or an active input-less circuit which is described by the following state space representation L : where z ∈ R n , A L ∈ R n×n , B L ∈ R n×1 , C L ∈ R 1×n . Clearly, L admits an input-output description by expressing the Laplace transform of v L as the product of the Laplace transform of i L and the equivalent impedance of L given by where I n is the identity matrix of order n and s is the complex variable.
The memristor MR is initially assumed to be an ideal flux-controlled memristor and thus described by the nonlinear flux-charge characteristic N : R → R, i.e.
where ϕ M and q M are the memristor flux and charge, whose standard definitions are recalled next [5]: In the voltage-current domain the memristor dynamics is modeled by the following nonlinear system where the derivative N (ϕ M ) is known as the memconductance of the memristor.
Since v L = v M and i L = −i M the dynamics of the circuit of Fig. 1 is described by the following equations C : Let z 0 (t) and ϕ 0 M (t) be the solution of (6) for t ≥ t 0 with initial conditions z 0 (t 0 ) = z 0 and ϕ 0 M (t) = ϕ M 0 , and let S C denote the set of all the solutions generated by varying z 0 ∈ R n and ϕ M 0 ∈ R. In the next section we determine under which conditions on L and MR, i.e., on the matrices A L , B L , C L and the memristor nonlinear characteristic N(·), there exists a one-to-one correspondence between S C and the solution set S of .

III. STRUCTURE OF THE DYNAMICALLY EQUIVALENT MEMRISTOR CIRCUIT
It has been shown that the state space of circuits containing an ideal memristors is composed of a continuum of invariant manifolds where either convergent or oscillatory and more complex dynamical behaviors can be displayed (see, e.g., [26] and references therein). These infinitely many invariant manifolds are parameterized by the manifold index I ∈ R whose value determines the specific invariant manifold where the dynamics is confined to lie. A general characterization of the invariant manifolds of the circuit of Fig. 1 C can be found in [37]. Here, we consider the case when A L is non-singular.
Proposition 3.1: Let A L be non-singular. Then, the invariant manifolds of the memristor circuit C are given by holds, i.e., the state space of the memristor circuit C is decomposed into a continuum of invariant manifolds.
Proof: Since A L is invertible, from the first equation of (6) we get and thus the second equation boils down tȯ ) is constant over time, thus showing that M I is an invariant manifold. To prove (8) it is enough to observe that for any fixed ϕ M all z ∈ R n belong to some manifold M I .

specific invariant manifold. Indeed, from the proof it follows the such a solution belongs to the invariant manifold M I with
Hence, the initial conditions z 0 , ϕ M 0 dictate the invariant manifold M I where the solutions are confined to lie. Proposition 3.1 states that the solution set S C of (6) can be obtained by collecting the sets of solutions belonging to M I for I ∈ R, i.e., We are now ready to determine the sought conditions ensuring that there is a one-to-one correspondence between S C and the solution set S in (2). The next result holds true.

Proposition 3.2: Let A be non-singular and suppose that the nonlinear characteristic of MR is chosen as
and the matrices of L are such that where T ∈ R n×n is a non-singular matrix. Then, for any x 0 ∈ R n and U 0 ∈ R the solution x 0 (t) of (1) for t ≥ t 0 with initial conditions x 0 (t 0 ) = x 0 and input u(t) = U 0 can be written as where z 0 (t) is the solution for t ≥ t 0 of (6) with initial conditions Moreover, z 0 (t) belongs to the invariant manifold M I with the index given by Proof: Consider the following equations which are obtained by first making the time-derivative of equations (1) and then replacingẋ and y with ξ and η, respectively. It can be readily verified that if x 0 (t), t ≥ t 0 , is the solution of (1) with initial conditions x 0 (t 0 ) = x 0 and input u(t) = U 0 , then the solution of (15) for t ≥ t 0 with initial conditions ξ(t 0 ) = Ax 0 + B f (Cx 0 ) + DU 0 and η(t 0 ) = Cx 0 is given by Now, the memristor circuit equations (6) with the nonlinear characteristic N(·) and the matrices A L , B L , C L are set according to (10) and (11), respectively, become By comparing (15) and (16) and rewriting the first equation of (16) as , it readily follows that the solutions z 0 (t) and ϕ 0 M (t) of (16) with the initial conditions in (13) are such that thus proving relation (12). The assumption that A is non-singular ensures that also A L is non-singular and thus from Proposition 3.1 we have that z 0 (t) belongs to an invariant manifold having the structure in (7). Taking into account (10) and (11), (7) can be rewritten as and hence, according to Remark 1, z 0 (t) belongs to the invariant manifold M I with By replacing z 0 with the expression in (13) and taking into account that ϕ M 0 = Cx 0 , we finally get Proposition 3.2 provides the conditions under which each solution x 0 (t) of is exactly reproduced via a suitable solution z 0 (t) of C , according to (12). We remark that, according to (10), the ideal memristor flux-charge characteristic needs to coincide with the function f (φ), which, as in the case of the FitzHugh model of Section IV, can be a monotonic odd function. When N(·) is non-monotonic and the memristor is active, it can be implemented in a standard way via a passive memristor in parallel to a negative resistor obtained using a current inverter (INIC). We also remark that an odd characteristic of the passive memristor may be obtained for instance by adequately biasing with some adapting circuitry PCM devices [18]. Another possibility is to resort to analog implementations via OTA or CMOS [13], [20] in combination with techniques to symmetrize the flux-charge characteristic as for instance using two memristors in a complementary resistive switching (CRS) configuration. The result in Proposition 3.2 can be strengthened as follows. Proposition 3.3: Let A be non-singular and let conditions (10)- (11) hold. If then there is a one-to-one correspondence between the solution sets S and S C of and C , respectively. Proof: If (17) is satisfied then (14) shows that there is a proportional relation between the index I of the invariant manifold and the value U 0 of the constant input. Hence, it is enough to prove that there exists a one-to-one correspondence between the solution sets S U 0 and S C I , with U 0 and I being related via (14). Recall that S U 0 is the set of all solutions of obtained for a fixed U 0 by varying the initial condition x 0 ∈ R n , while S C I is the set of all solutions of C with initial conditions z 0 , ϕ M 0 belonging to the following set Since L I ≡ R n , to prove one-to-one correspondence between S and S C it is enough to show that there not exist two different solutions in S to which correspond the same solution in S C .
have the initial conditionsz 0 ,φ M 0 andẑ 0 , ϕ M 0 , respectively. According to (13), we havē if and only ifz 0 =ẑ 0 andφ M 0 =φ M 0 . These two conditions can be equivalently written as and respectively. Exploiting (19), condition (18) boils down to and thus it is solved only ifx 0 =x 0 , which is a contradiction. Remark 2: The proof of Proposition 3.3 shows that each invariant manifold M I replicates all the solutions of (1) obtained for a fixed constant input u(t) = U 0 once I and U 0 satisfy the linear relation (14). In the case that displays different attractors for different values of U 0 , this relation permits to keep track of the invariant manifolds where these attractors are reproduced. In [37] it has been shown that the introduction of pulse programmed sources in the memristor circuit makes it possible to switch the circuit dynamics between different attractors. This suggests that the case of non-constant input u can be tackled by incorporating suitably designed voltage/current sources in the two-terminal element L of the circuit of Fig. 1. This issue will be investigated in Section V.
Proposition 3.2 defines the structure that equations (3) and (5) should enjoy to reproduce the solutions of system . Specifically, the flux-controlled memristor is completely designed once the nonlinear characteristic satisfies (10), while different matrices A L , B L , and C L can be chosen by varying the matrix T .
This degree of freedom is used to select a specific circuital realization of the equations (3) governing the two-terminal element. We first observe that for any choice of A L , B L , and C L the impedance (4) of L satisfies Hence, we have to synthesize either a passive or an active two-terminal element L having such an impedance. This can be done by exploiting one of the available procedures (see, e.g., [42] and [43]), together with the impedance scaling technique. Note that Assumption 1 ensures that the strictly proper rational function L(s) has order equal to n. This implies that L can be synthesized with n C ≥ 0 capacitors and n L ≥ 0 inductors with n = n C + n L . Let v C i , i = 1, . . . , n C , and i L j , j = 1, . . . , n L , denote the voltages of the capacitors and the currents of the inductors, respectively, and consider the by applying the Kirchhoff's law to the loops and nodes of the synthesized L, we can write the following state space representation for some suitable matricesĀ L ,B L , andC L . Since by con-structionC there exists a non-singular matrixT ∈ R n×n satisfying (11), i.e., solving the linear equations Clearly, the designed circuit obeys equations (6) with A L = A L , B L =B L , C L =C L , and N(·) = − f (·) and hence, according to Proposition 3.2, it is capable to reproduce all the solutions of . Specifically, the solution x 0 (t) of (1) for t ≥ t 0 with initial conditions x 0 (t 0 ) = x 0 and input u(t) = U 0 is given by (12) once T is replaced byT and z 0 (t) is the solution generated by the designed circuit for t ≥ t 0 with initial conditions as in (13) with T −1 replaced byT −1 .

A. The Case of Singular Matrix A
In this subsection we show how the assumption of A being non-singular can be relaxed. Indeed, if the matrix A of is non-invertible we can rewrite equations (1) as where k ∈ R is a parameter to be chosen such that det(A−k BC) = 0. Clearly, for such a k all the results obtained in the case of a non-singular A are still valid once A and f (y) are replaced by A k . = A − k BC and f k (y) . = f (y) + ky. In particular, it can be checked that the impedance of the two-terminal element L to be synthesized becomes where L(s) is as in (20). This relation between L k (s) and L(s) is exploited to show that it is possible to find k such that A k is non-singular. The next result holds true. and L k (s) can be written as Memristor circuit capable to reproduce the dynamics of the FitzHugh-Nagumo neuron model: the passive two-terminal elementi L (blue box) with R m , L m , and C m as in (32), and the flux-controlled memristor (red box) with the nonlinear characteristic (30).
where the second equality follows from (25). Hence, we have the following equality

IV. APPLICATION TO THE FITZHUGH-NAGUMO MODEL
In this section the application of the previous results will be illustrated by designing a memristor circuit capable to reproduce the dynamics of the FitzHugh-Nagumo neuron model [44], [45], [46], [47], which is a second order approximation of the celebrated Hodgkin-Huxley model [48]. Some other memristive FitzHugh-Nagumo circuits have been already proposed for different purposes [49], [50], [51]. The model FitzHugh-Nagumo equations read where V represents the membrane potential, W is an internal variable, the coefficients a, b, c are positive tuning parameters, and I is a constant forcing current, used to stimulate the neuron response. It can be readily verified that model (27) admits the state space representation in (1) once Since a, b, c > 0, we have and hence A is non-singular and condition (17) holds.
To design the memristor circuit of Fig. 1, we first observe that the flux-controlled memristor MR should be implemented according to equation (5) where the nonlinear characteristic must satisfy (10), which implies Exploiting the expression of A, B, and C in (29), we get that the impedance (4) of the two-terminal element L is given by L(s) = s + ac s 2 + acs + ab (31) and, therefore, the corresponding admittance is It can readily verified that this admittance can be synthesized via the parallel interconnection of a capacitor C m with the series of an inductor L m and a resistor R m , once Figure 2 depicts the synthesized circuit whose governing equations are ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩v where To apply Propositions 3.2 and 3.3 it remains to compute the matrixT solving conditions (23). Indeed, we get where T i j , i = 1, 2, j = 1, 2, are the entries ofT . The above equations boil down to which, taking into account (32), are solved for According to Proposition 3.3, there is a one-on-one correspondence between the solutions of (27) generated by the constant forcing current I and those of (33) belonging to the invariant manifold with the manifold index given by Moreover, from Proposition 3.2 and in particular from (13), it follows that the solution of (27) with initial condi-

A. Numerical Simulations
Consider the FitzHugh-Nagumo model (27). Without any loss of generality, both the state variables and the time variable can be assumed to be expressed in normalized values, which are obtained from the corresponding biological values by shifting the axes origin and applying a scale factor. Then, consider the normalized parametric configuration   which is coherent with the original system (see [44], [45], and [46]). The system is numerically simulated in the Matlab environment using the ode45 solver for different constant values of the input I and of the initial conditions x 0 . Figure Figure 4 illustrates the same trajectories in their temporal evolution.
Consider now the corresponding memristor circuit depicted in Fig. 2. From (32), the normalized equations of circuit (33) feature In view of practical implementation, an impedance scaling technique has to be employed to obtain realistic values of the circuit components. The previous starting points x (1) 0 , x  and (z (2) 0 , ϕ (2) M 0 ) belong to the same manifold of index I = 0.16, while (z (3) 0 , ϕ (3) M 0 ) belongs to that of index I = −0.64. A picture of the corresponding manifolds is provided by Fig. 5. The trajectories of the circuit originated from those starting points are shown in Fig. 6.

V. EXTENSION TO NON-CONSTANT INPUTS u
In this section we relax the assumption that the input u of is constant. We recall that typically displays an attractor for many constant values of the input u, as illustrated in Fig. 7 for the FitzHugh-Nagumo model. Consider the case where exhibits two different attractors for the input values U (1) 0 and U (2) 0 , U (1) 0 = U (2) 0 , and let u(t) be varied between these values in order to make the dynamics switching between the different attractors.
As an example, we can consider the input where r (t) satisfies r (t 1 ) = U (1) 0 and r (t 2 ) = U (2) 0 . With some abuse of notation, we denote by x 0 (t) the solution with initial condition x(t 0 ) = x 0 and input u as in (38). We suppose that at t = t 1 such a solution lies onto the attractor corresponding to the constant input U (1) 0 , while at t = t 2 it belongs to the stability region of the attractor relative to the constant input U (2) 0 so that it converges towards this attractor for t > t 2 . Propositions 3.2 and 3.3 ensure that the attractors Fig. 6. State space orbits of the FitzHugh-Nagumo circuit (33)-(37) originated from the initial conditions x (1) 0 , x (2) 0 , x (3) 0 transformed according to (35). The circle, triangle and square correspond to the equivalent symbols used in Fig. 3 for the simulation of the FitzHugh-Nagumo system. The trajectories starting from the circle and triangle lie on the manifold of index I = 0.16, while the orbit originated from the square lies onto the manifold of index I = −0.64.  (2) 0 . Clearly, the switching solution x 0 (t) cannot be reproduced by C since its dynamics cannot move from one invariant manifold to another. However, it can be shown that x 0 (t) can be replicated by introducing suitable voltage/current sources in the two-terminal element L which has been synthesized in the case of constant inputs. We recall that in this case the state space representation and the impedance are as in (21) and (22), respectively, and that T =T in Proposition 3.2. Proceeding as in [37], we replace each capacitor C i with the same capacitor with in parallel a current source I (s) i , i = 1, . . . , n C , and each inductor L j with the same inductor with in series a voltage source V (s) j , j = 1, . . . , n L . It follows that the state space representation of the so-modified L becomes where ∈ R n is the vector of the voltage and current sources, i.e., and ∈ R n×n is the following non-singular diagonal matrix Hence, by joining equations (39) with (5) and exploiting conditions (10) and (11), we get that the dynamics of the modified memristor circuit obeys The next result is instrumental for the design of the vector ensuring that equations (42) replicate the switching solution of . Proposition 5.1: Let x 0 (t) denote the solution of with initial condition x 0 (t 0 ) = x 0 and non-constant input u(t) such thatu(t) is piecewise-continuous for t ≥ t 0 , and letz 0 (t) and ϕ 0 M (t) be the solutions for t ≥ t 0 of (42) with initial conditions Suppose that the vector of the current and voltage sources satisfies for t ≥ t 0 . Then, x 0 (t) can be written as Proof: By proceeding as in the proof of Proposition III we have that the equations admit the solutions ξ(t) =ẋ 0 (t) and η(t) = Cx 0 (t) for the initial condition ξ(t 0 ) = Ax 0 + B f (Cx 0 )+ Du(t 0 ) and η(t 0 ) = Cx 0 . We observe that equations (42) can be rewritten as once is replaced with its expression in (44). By comparing (47) and (46), it readily follows that the solutions z 0 (t) and ϕ 0 M (t) of (47) with the initial conditions in (43) are such that ϕ 0 M (t) = Cx 0 (t) andTz 0 (t) =ẋ 0 (t), thus proving (45). (38) then u(t 0 ) = U (1) 0 , u(t 2 ) = U (2) 0 and (t) = 0 for both t ∈ [t 0 , t 1 ] and t > t 2 . According to Proposition III, this implies that the solutionsz 0 (t) and ϕ 0 M (t) of (42) is confined to lie onto the invariant manifold

Remark 4: It is worth noting that if u(t) is as in
is non-constant,z 0 (t) and ϕ 0 M (t) are not belonging to any invariant manifold. This means that the manifold index is no longer constant. In fact, we havė where the second equality follows by replacingTż 0 (t) anḋ ϕ 0 M (t) with the expressions in the right-side terms of (47), respectively. Since (48) we get Hence, for t ∈ [t 1 , t 2 ] the index of the manifold varies linearly with the non-constant input of .

A. The FitzHugh-Nagumo Model
According to the above results, the implementation of the memristor circuit when input u is non-constant is achieved by introducing as many voltage and current sources as the number of inductors and capacitors, respectively, in the circuit implementation already obtained for the constant input case. Taking into account that all the voltage sources must be in series of the inductors, and the current sources in parallel of the capacitors, the application of this procedure to the circuit equations (33) brings to the circuit implementation illustrated   in Fig. 8, which is described by the following system: Observe that the vector of the sources is = (I (s) 1 , V (s) 1 ) , just in agreement with (40). Also, note that = diag (C m , L m ) in accordance to (41). In order to change the manifold index as in (49), condition (44) must be enforced. Taking into account (32), it turns out that Therefore, the voltage source V (s) 1 is not necessary to steer the manifold index. From (49) it follows that such a enforces in the manifold index the following dynamics: where u stands for the time-varying stimulus of the FitzHugh-Nagumo system, as described by (28)- (29). Note that its derivativeu is the real objective of the design, since it defines the actual implementation of the current source, that, in this example, has been chosen as follows: In Fig. 9 the circuit (50) is numerically simulated, as before, by the Matlab software, and the same configuration (37)  1 source (or equivalentlẏ u) is used to steer the trajectory from a manifold to another. In particular, it is constant and non-negative in (95, 105) and it acts as a linear ramp in (240, 260). Therefore, the manifold index varies as linear ramp in the first period, and as a decreasing parabolic ramp in the second one. Out of these intervals the source is constant at zero, and therefore the manifold index remains unchanged. Figure 10 illustrates the same resulting trajectory in the state space of the circuit, while the dynamics of the manifold index is reported in Fig. 11.

B. The Duffing System
The application of the previous results will be illustrated by designing a memristor circuit capable to reproduce the dynamics of the forced Duffing system. Other memristive Duffing type circuits have been proposed in the recent literature (see, e.g., [52] and [53]), but their design has been mainly inspired  [43] with R, C 1 , and C 2 as in (55), and the flux-controlled memristor (red box) with the nonlinear characteristic (54).
0 , and x (4) 0 (the circle symbols) when the system if respectively forced by u 1 (t), u 2 (t), u 3 (t), and u 4 (t). Second row: Different dynamics exhibited by the Duffing circuit of Fig. 12 for configurations equivalent to those reported in the left column for the Duffing system. E, F, G, H: State space orbits of originated from the initial conditions (57) computed on the couples (x (1) 0 , u 1 (t 0 )), (x (2) 0 , u 2 (t 0 )), (x (3) 0 , u 3 (t 0 )), and (x (4) 0 , u 4 (t 0 )) derived from the original configurations of the Duffing system. In each case the circuit is driven by the related I  by the original nonlinear circuit [54]. The Duffing system equation reads [55] ξ(t) + δξ (t) + αξ + βξ 3 where the forcing input u(t) is any non-constant input satisfying the assumptions of Proposition 5.1. System (51) can be put in form (1) by defining and the corresponding L(s) is For the implementation's sake and without any loss of generality, consider the equivalent though extended form (24), where k plays the role of an additional degree of freedom. Then, the synthesis of the circuit can be performed starting from impedance (25), which turns out equal to L k (s) = β Indeed, if α + kβ > 0, the impedance (53) can be obtained via the circuit L depicted in Fig. 12, where each capacitor has already been coupled with as many parallel current sources. The equations governing the two-terminal element L are ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩v and, thus, they admit the form (39) oncez = (v C 1 , v C 2 ) and The complete circuit is obtained by introducing the dynamics (5) of the memristor MR with The solution of (23) leads to and also to the inverse relationships between the circuit elements and the original dynamics' parameters, which are Finally, observe that equation (44) gives and then the current source I (s) 1 turns out unnecessary. Indeed, the manifold index is related to the input u as follows: The Duffing system (51) in its state representation (52) has been numerically simulated in the Matlab environment for different inputs u(t) when the parametric configuration is α = −1, β = 1, δ = 0.3. Figures 13.A, 13.B, and 13.C shows the trajectories generated from x Memristor circuit capable to reproduce the dynamics of the FitzHugh-Nagumo neuron model with a charge-controlled memristor (red box) of nonlinear characteristic (60) and a passive two-terminal elementi L (blue box) with R m , L m , and C m as in (61).  (1.2t), respectively. As the input amplitude increases, Duffing system exhibits a limit cycle undergoing a period doubling cascade, which ends into chaos. The chaotic attractor for u 4 (t) = 0.50 cos(1.2t) is reported in Fig. 13.D for the initial condition x (4) 0 = (−0.8553, 0.5371) . According to (55) and assuming k = 2 so that α + kβ = 1, the corresponding memristor circuit features the parameters R = 0.5, C 1 = 0.3, C 2 = 13.3333. The starting points of the simulations of Fig.s 13.A, 13.B, 13.C, and 13.D are transformed into initial conditions of the circuit according to (43), which reads

VI. EXTENSION TO CHARGE-CONTROLLED MEMRISTORS
In this section we consider the case when the memristor MR in Fig. 1 is assumed to be an ideal charge-controlled memristor. Such a memristor is described by a nonlinear charge-flux characteristic which is denoted by N(·) as in the case of the flux-controlled memristor. It relates the memristor charge and flux as follows In the voltage-current domain the memristor dynamics is modeled by the following nonlinear system where the derivative N (q M ) is known as the memristance of the memristor. Note that in the case of charge-controlled memristor i M and v M are the input and the output of MR, respectively. Hence, to ensure that the memristor circuit of Fig. 1 is well-posed, the two-terminal element L should have v L as input and i L as output, i.e., L is either a passive or an active input-less circuit which is described by the following state space representation L : where z ∈ R n , A L ∈ R n×n , B L ∈ R n×1 , C L ∈ R 1×n . Observe that in this case the strictly proper rational function is the admittance of the two-terminal element.
Since v L = v M and i L = −i M , it follows that the dynamics of the circuit of Fig. 1 obeys the following equations C : which have been labeled with the same symbol C used for the equations (6) in the flux-controlled memristor case. By comparing (6) and (59) we note that the unique difference is that ϕ M has been replaced with q M . Hence, we can conclude that all the results obtained for the flux-controlled memristor are still valid for the case of the charge-controlled memristor, once the flux ϕ M is replaced by the charge q M . Clearly, now L should be synthesized by ensuring that its admittance is equal to (4).

A. The Charge-Controlled FitzHugh-Nagumo Circuit
In this subsection a charge-controlled memristor circuit is derived for the FitzHugh-Nagumo neuron model by following the previous procedure. In this case, the memristor is described by (58) along with the nonlinear characteristic while L(s) in (31) represents the admittance of the two-terminal element L between v L and i L . This implies that the corresponding impedance can be written as L −1 (s) = s 2 + acs + ab s + ac = s + ab s + ac .
Hence, it turns out that it can be synthesized by the series of an inductor L m with the parallel connection of a resistor R m and a capacitor C m , once Figure 14 depicts the implemented circuit which is governed by the following equations ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩v where v M =Ṅ (q M ) and i M = −i L m . Since v L = v M and i L = −i L m , we get that (62) can be written as in (6) once z = (v C m , i L m ) and A L =Ā L , B L =B L , C L =C L with Finally, the matrixT solving conditions (23) is

VII. CONCLUSION
In this paper the dynamics of a class of forced nonlinear systems which includes several systems displaying complex dynamics, has been considered. First, it has been shown that there exists a one-to-one correspondence between the dynamics of the nonlinear system obtained for a given constant value of the forcing input and that displayed on one of the invariant manifolds of a suitably synthesized input-less circuit composed of a two-terminal element connected with an ideal flux-controlled memristor. In particular, it turns out that in the case of the FitzHugh-Nagumo model there is a linear relation between the values of the injected current and the index of the invariant manifolds. Then, it is shown that, even in the case of a non-constant forcing input, the system dynamics can be still replicated by introducing suitably programmed voltage/current sources in the two-terminal element of the circuit, as illustrated via the application to the FitzHugh-Nagumo model and the Duffing system. Also, it is observed that quite similar results hold if the flux-controlled memristor is replaced by a chargecontrolled one. Several future research issues can be foreseen, such as the practical implementation of ideal memristor and their interconnection with suitably designed two-terminal elements, as well as the study of the capability of circuits with real physical memristors to replicate the dynamics of known oscillatory systems.