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.

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.

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 inmemory computing, i.e., the co-location of computation and memory in the same device (Ielmini and Wong, 2018). Inmemory 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., 2016Ascoli 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 fluxcharge analysis method (FCAM), has been introduced for the analysis of memristor circuits in the flux-charge domain (Corinto and Forti, 2016. 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, 2000Izhikevich, , 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. 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.

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.
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φ : R → R of the memristor relating the charge q M and the flux ϕ M is assumed to have both a linear and a cubic term. Specifically, where the constant terms s 0 and s 1 satisfy Throughout the paper, we assume that all the circuit parameters R, L, C, s 0 , s 1 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. It can be readily verified that the circuit dynamics is described by the following third-order system of differential equations where D is the time-derivative operator 1 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 Hence, in the input-less case the state space of the memristor circuit is decomposed into infinitely many invariant manifolds of the form where Q 0 is the index of the manifold whose value depends on the circuit initial conditions according to the relation Q 0 = q M (t 0 ) + Cv C (t 0 ). Note that the invariant manifolds M(Q 0 ) are planar surfaces parameterized by the manifold index Q 0 , as illustrated in Figure 3.
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 0 ) can be readily singled out. Since the first equation of Equation (3) has been used to obtain M(Q 0 ), it follows that the dynamics is characterized by the two Hence, the dynamics onto M(Q 0 ) obeys the second-order system of differential equations 1 Throughout the paper, Df (t) denotes the derivative of the function f (t) at t =t.

A B
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.
FIGURE 3 | Invariant manifolds of the input-less circuit: the state space trajectory generated by the solution of Equation (3) 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 0 ). The dynamical analysis of Equation (7) for any value of the index manifold Q 0 will be pursued in section 3. Note that Q 0 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 0 , and s 1 .

INPUT-LESS MEMRISTOR CIRCUIT DYNAMICS
In this section, we summarize the properties of the dynamics onto M(Q 0 ) 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.

Equilibrium Points
It readily follows that each invariant manifold M(Q 0 ) has a unique equilibrium point at (q M , i L ) = (Q 0 , 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 , 0) for any value of the index Q 0 . 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 0 ), we compute the Jacobian J(Q 0 ) of (7) at (q M , i L ) = (Q 0 , 0), getting The eigenvalues of J(Q 0 ) have a negative real part if |Q 0 | > Q 0 and a positive real part if |Q 0 | <Q 0 , while they are pure imaginary at Q 0 = ±Q 0 . This implies that the equilibrium point at (q M , i L ) = (Q 0 , 0) is ensured to be (locally) asymptotically stable if |Q 0 | >Q 0 and unstable if |Q 0 | <Q 0 . Hence, the equilibrium point of each manifold M(Q 0 ) with |Q 0 | >Q 0 is an attractor of the input-less memristor circuit. Figure 4 summarizes the dynamical behavior around the equilibrium point of M(Q 0 ).

Limit Cycles
To investigate the presence of limit cycles on M(Q 0 ) we resort to the DF method, a classical technique within the HB FIGURE 4 | Stable (green half lines) and unstable (red segment) equilibrium points: the trajectories starting on the planes with Q 0 < −Q 0 and Q 0 >Q 0 converge toward the corresponding equilibrium points; the trajectory starting on the plane with Q 0 ∈ (−Q 0 ,Q 0 ) converges toward the stable limit cycle (green).
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 inputoutput relation where L(D) and L 0 (D) are suitable rational functions.
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 = Dq M the second equation of Equation (7) becomes Taking into account that −s 0 Dq M (t)+s 1 q 2 M (t)Dq M (t) = Dφ M (t), the above differential equation can be rearranged in the following equivalent input-output form Hence, system (7) can be equivalently described via the inputoutput relation (10) with L(D) and L 0 (D) given by 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., The corresponding periodic output ϕ M (t) of period 2π/ω of the non-linear subsystem can be written as where N 0 (A, B) is the bias gain and N 1 (A, B) is first harmonic gain, which are known as describing function terms, while ϕ h (t) contains the higher order harmonics, i.e., with γ h and θ h , h = 2, . . ., being suitable constants. The final step of the DF method consists in first computing the Frontiers in Neuroscience | www.frontiersin.org periodic signal of period 2π/ω given by the sum of the outputs y(t) and y 0 (t) of the two linear subsystems driven by −ϕ M (t) and the constant signal Q 0 , 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)Ae 0t and N 1 (A, B)B(e jωt + e −jωt )/2, respectively, the periodic signal given by the sum of y(t) and y 0 (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 2 . The so computed periodic signal has the following expression (17) 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 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 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 (0) = L(D) = b n−1 D n−1 + . . . + b 1 D + b 0 D n + a n−1 D n−1 + . . . + 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 ∈ R, σ ∈ 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. (1) and (15) it can be verified that the gains N 0 (A, B) and N 1 (A, B) have the following expressions:

1, while from Equations
Then, the Equation (18) reduces to and the complex equation (19) boils down to which can equivalently be written by separating the real and imaginary parts as Equations (22), (24), and (25) are then solved for A =Â = Q 0 , Hence, for each index Q 0 such that |Q 0 | <Q 0 there exists a PLC of the following form (27) This implies that each manifold M(Q 0 ) contains a unique PLC for |Q 0 | <Q 0 . Since the trajectory described by the PLC in the (q M , i L )-plane lies on the following ellipse Note that the ellipse is centered at the equilibrium point (Q 0 , 0) and its size is maximum at Q 0 = 0 and decreases to zero as the index Q 0 tends toQ 0 , thus collapsing to the equilibrium point. Moreover, since for |Q 0 | <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 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 0 . The red circles denote the unstable equilibrium points.
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 0 such that |Q 0 | >Q 0 and a stable PLC for each value of the index Q 0 such that |Q 0 | <Q 0 . Figure 6 illustrates this multistability scenario for the following normalized values of the circuit parameters R = 0.4 , C = 0.1 , L = 1.5 , s 0 = 0.7 , s 1 = 0.3 , (30) from which it follows thatQ 0 = 1.
In particular, it turns out that at Q 0 = ±Q 0 = ±1 a typical behavior of (supercritical) Hopf bifurcation is observed. However, since it is obtained by varying the index Q 0 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 0 | < 1 the corresponding invariant manifold M(Q 0 ) 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 = 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 0 .
It can be readily checked that the minimum and the maximum of the first-order approximation are given by respectively. We finally observe that similar results can be derived also for values of the circuit parameters different from those in Equation (30).

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 0 ) contains either a unique stable equilibrium point or a stable limit cycle. Moreover, an approximation of the limit cycle belonging to M(Q 0 ), with Q 0 such that |Q 0 | <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 0 , 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 It turns out that (3) reduces to the third-order system where Q, q M , and i L are the new state variables. Suppose that at time t = t i ≥ t 0 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 0 , t i ] and t ∈ [t i + , +∞) and such that From the first equation of Equation (34) it follows that which implies Q(t) = Q (i) 0 for t ∈ [t 0 , t i ] and Q(t) = Q (i) 0 + for t ∈ [t i + , +∞). For instance, in the case of a rectangular pulse we have Q(t) = Q (i) 0 + (A/ )(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 M(Q (i) 0 ) for t ∈ [t 0 , t i ], it moves toward M(Q (i) 0 + ) during the interval [t i , t i + ], it reaches M(Q (i) 0 + ) 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 M(Q (i) 0 ) has a stable equilibrium point and in particular the index Q (i) 0 is negative and hence Q (i) 0 < −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 (i) 0 , Q (i) 0 , 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 M(Q (i) 0 + ) still contains a stable equilibrium point if < |Q (i) 0 | −Q 0 , while it contains a stable limit cycle if |Q (i) 0 | −Q 0 < < |Q (i) 0 | +Q 0 and again a stable equilibrium point if > |Q (i) 0 | +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).
Specifically, after reaching M(Q (i) 0 + ) at time t = t i + at the point (Q(t i + ), q M (t i + ), i L (t i + )) with Q(t i + ) = Q (i) 0 + , 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 M(Q (i) 0 + ), 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 (i) 0 , 0) and Q(t) = Q (i) 0 + (t −t i )/ . 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 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 M(Q (i) 0 ) to the manifold M(Q (i) 0 + 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.
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 M(Q (i) 0 ); at time t = t i + the dynamics reaches M(Q (i) 0 + ) and the memristor charge q M starts displaying a behavior similar to that reported in Figure 2B; at time t = t i + A B 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)  + T an inverse current pulse is generated in order to make the dynamics come back to the stable equilibrium point in M(Q (i) 0 ) (reset mechanism).
Hence, the problem is how the parameters Q (i) 0 , , , 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 M(Q (i) 0 + ) we have that q M and i L are quite close to Q (i) 0 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 M(Q (i) 0 + ) is obtained by putting Q 0 = Q (i) 0 + 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 + ), i L (t i + )) ≈ (Q (i) , 0), it is possible to set (q M (t i + ), i L (t i + )) as the point of the manifold M(Q (i) 0 + ) 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 M(Q (i) 0 ). Hence, the resulting behavior of the memristor charge q 0 M (t) can be expressed as and is depicted in Figure 10. Note that the parameter T has been chosen as the period of the PLC, i.e.: Hereafter, we denote by i 0 L (t) the time-derivative of q 0 M (t), i.e.: Now, to set (q M (t i + ), i L (t i + )) as the point of the manifold M(Q (i) 0 + ) where the PLC has its minimum value (31), the following condition must be satisfied with Q 0 = Q (i) 0 + . It is not difficult to verify that there are many solutions Q (i) 0 and Q 0 of (41) such that Q (i) 0 < −Q 0 and |Q 0 | <Q 0 . We are interested in the one with the minimum value of Q (i) 0 because this choice ensures that the stable equilibrium point on the initial manifold M(Q (i) 0 ) is at the maximum distance from the invariant manifold M(Q 0 ) where the Hopf bifurcation happens, thus guaranteeing a larger robustness margin against spurious pulses of small area. The minimum value of Q (i) 0 can be and since from Equations (41) and (9) we have Moreover, since Q 0 in Equation (41) stands for Q (i) 0 + A, from Equation (42) we finally get 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 where σ lies in the interval [0.01, 0.05]. Summing up, to obtain q 0 M (t) in Equation (38), the parameters T, Q (i) 0 , , 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 0 , and s 1 . Hence, by suitably adapting these parameters it is possible to vary the features of q 0 M (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 (i) 0 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 (i) 0 , 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) and choosing σ = 0.02 in Equation (45) we obtain = 0.0487 .
The corresponding predicted behavior q 0 M (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 0 M (t), i 0 L (t))) are depicted in Figure 11B.
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 0 M (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.

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.