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.

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 φ^:ℝ→ℝ 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₀ and s₁ satisfy

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.

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¹ 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

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.,

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₀ 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.

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 1C(Q0-qM(t)). Hence, the dynamics onto M(Q₀) obeys the second-order system of differential equations

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₁.

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.

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 |Q0|<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

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₀ 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

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

which implies Q(t)=Q0(i) for t ∈ [t₀, t_i] and Q(t)=Q0(i)+Λ for t ∈ [t_i + Δ, +∞). For instance, in the case of a rectangular pulse we have Q(t)=Q0(i)+(A/Δ)(t-ti).

By comparing the second and third equations of Equation (34) with those of Equation (7), it follows that the circuit dynamics lies onto M(Q0(i)) for t ∈ [t₀, t_i], it moves toward M(Q0(i)+Λ) during the interval [t_i, t_i + Δ], it reaches M(Q0(i)+Λ) 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(Q0(i)) has a stable equilibrium point and in particular the index Q0(i) is negative and hence Q0(i)<-Q-0. Moreover, we assume that at t = t_i the circuit state is at the equilibrium point, which means that (Q(ti),qM(ti),iL(ti))=(Q0(i),Q0(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 M(Q0(i)+Λ) still contains a stable equilibrium point if Λ<|Q0(i)|-Q-0, while it contains a stable limit cycle if |Q0(i)|-Q-0<Λ<|Q0(i)|+Q-0 and again a stable equilibrium point if Λ>|Q0(i)|+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(Q0(i)+Λ) at time t = t_i + Δ at the point (Q(t_i + Δ), q_M(t_i + Δ), i_L(t_i + Δ)) with Q(ti+Δ)=Q0(i)+Λ, 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(Q0(i)+Λ), 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 (qM(ti),iL(ti))=(Q0(i),0) and Q(t)=Q0(i)+Λ(t-ti)/Δ. 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(Q0(i)) to the manifold M(Q0(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.

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(Q0(i)); at time t = t_i + Δ the dynamics reaches M(Q0(i)+Λ) 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 M(Q0(i)) (reset mechanism).

Hence, the problem is how the parameters Q0(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 M(Q0(i)+Λ) we have that q_M and i_L are quite close to Q0(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 M(Q0(i)+Λ) is obtained by putting Q0=Q0(i)+Λ 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 (qM(ti+Δ),iL(ti+Δ))≈(Q(i),0), it is possible to set (q_M(t_i + Δ), i_L(t_i + Δ)) as the point of the manifold M(Q0(i)+Λ) 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(Q0(i)). Hence, the resulting behavior of the memristor charge qM0(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 iL0(t) the time-derivative of qM0(t), i.e.:

Now, to set (q_M(t_i + Δ), i_L(t_i+Δ)) as the point of the manifold M(Q0(i)+Λ) where the PLC has its minimum value (31), the following condition must be satisfied

with Q0=Q0(i)+Λ. It is not difficult to verify that there are many solutions Q(i)0 and Q₀ of (41) such that Q0(i)<-Q-0 and |Q0|<Q-0. We are interested in the one with the minimum value of Q0(i) because this choice ensures that the stable equilibrium point on the initial manifold M(Q0(i)) 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 Q0(i) can be computed in closed form by minimizing the right hand side of Equation (41). It turns out that

and since

from Equations (41) and (9) we have

Moreover, since Q₀ in Equation (41) stands for Q0(i)+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 qM0(t) in Equation (38), the parameters T, Q0(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 qM0(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 Q0(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, Q0(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

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

The corresponding predicted behavior qM0(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),qM0(t),iL0(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 qM0(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.

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.

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.

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

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.