Electromechanical memcapacitive neurons for energy-efficient spiking
  neural networks
In this article, we introduce a new nanoscale electromechanical device -- a
leaky memcapacitor -- and show that it may be useful for the hardware
implementation of spiking neurons. The leaky memcapacitor is a movable-plate
capacitor that becomes quite conductive when the plates come close to each
other. The equivalent circuit of the leaky memcapacitor involves a
memcapacitive and memristive system connected in parallel. In the leaky
memcapacitor, the resistance and capacitance depend on the same internal state
variable, which is the displacement of the movable plate. We have performed a
comprehensive analysis showing that several spiking types observed in
biological neurons can be implemented with the leaky memcapacitor. Significant
attention is paid to the dynamic properties of the model. As in leaky
memcapacitors the capacitive and leaking resistive functionalities are
implemented naturally within the same device structure, their use will simplify
the creation of spiking neural networks.
Zixi Zhang, Yuriy V. Pershin, Ivar Martin
# Electromechanical memcapacitive neurons for energy-efficient spiking neural networks

###### Abstract

In this article, we introduce a new nanoscale electromechanical device - a leaky memcapacitor - and show that it may be useful for the hardware implementation of spiking neurons. The leaky memcapacitor is a movable-plate capacitor that becomes quite conductive when the plates come close to each other. The equivalent circuit of the leaky memcapacitor involves a memcapacitive and memristive system connected in parallel. In the leaky memcapacitor, the resistance and capacitance depend on the same internal state variable, which is the displacement of the movable plate. We have performed a comprehensive analysis showing that several spiking types observed in biological neurons can be implemented with the leaky memcapacitor. Significant attention is paid to the dynamic properties of the model. As in leaky memcapacitors the capacitive and leaking resistive functionalities are implemented naturally within the same device structure, their use will simplify the creation of spiking neural networks.

## I Introduction

Information processing in biological systems relies on a complex network of interacting neurons. Each neuron, when subjected to a stimulus responds by outputting a signal that typically has a form of unharmonic spikes. A variety of models of the spiking phenomenon have been proposed, most famously the Hodgkin-Huxley (HH) model, which successfully captured many observed features of spiking in neuronal membranes [14; 17; 18; 19; 20; 21; 24]. This model prompted the development of bio-chemically based information processing models in the subsequent decades [24; 5; 40; 27].

HH model attempts to accurately capture the properties of ion channels with memory and consequently is quite complex. A number of simplified models, like the integrate-and-fire model, appeared that can achieve a variety of spiking behaviors in response to different types of stimulation [14; 12; 14; 23; 24; 28; 4]. Despite being distinct from HH, one may expect that the qualitative features such as instabilities/bifurcations, limit cycles, and synchronization among neurons are quite universal and robust due to the general underlying principles of the dynamical systems [39]. The simplified models, however, may have the important advantage of being easier to implement artificially using existing materials and devices. This leads to an exciting possibility of biologically-inspired _neuromorphic_ information processing systems implemented in the solid-state.

In the above context, the class of memory circuit elements [11] becomes increasingly important because of the capacity of memory circuit elements to store and process the information on the same physical platform [10]. The memory circuit elements (in a pure form) are resistors, capacitors, and inductors with memory whose response is defined by the equations

\[y(t) =g\left(\mathbf{x},u\right)u(t)\ \, \tag{1}\] \[\dot{\mathbf{x}} =\mathbf{f}\left(\mathbf{x},u\right)\ \, \tag{2}\]

where \(y(t)\) and \(u(t)\) are any two complementary circuit variables (i.e., current, charge, voltage, or flux), \(g\left(\mathbf{x},u\right)\) is a generalized response, \(\mathbf{x}\) is a set of \(n\) state variables describing the internal state of the device, and \(\mathbf{f}\left(\mathbf{x},u\right)\) is a continuous \(n\)-dimensional vector function. Depending on the choice of the complementary circuit variables, Eqs. (1) and (2) are used to define memristive, memcapacitive, or meminductive systems [11]. Examples of physical realizations of memory circuit elements can be found in the review paper [35].

For more than a decade, significant attention has been paid to the application of memristive systems in neuromorphic computing. Indeed, the memristive systems share several common characteristics with biological synapses such as
the two-terminal structure, adaptivity, and high integration density. One of the first works in this area was the demonstration of associative memory with memristive neural networks by one of us [34]. Interestingly, the relevance of the memristive equations to the HH model was established already in 1976 by Chua and Kang [7; 8]. Specifically, they proposed that "the potassium channel of the Hodgkin-Huxley model should be identified as a first-order time-invariant voltage-controlled memristive one-port and the sodium channel should be identified as a second-order time-invariant voltage-controlled memristive one-port" [8]. Even though most of the existing research focuses on purely electronic schemes, there are also models that rely on the mechanical realization of memory and spikes [6; 13; 15; 22; 25].

Neuromorphic applications of memcapacitive systems have received much less attention. This may be explained by the fact that in general memcapacitive devices [30; 31; 32; 33] have been much less studied (in comparison to the memristive ones). At the same time, the reactive nature of the capacitive response is very promising for low-power computing applications as at equilibrium the memcapacitive devices do not consume any power. Thus, currently, the application of memcapacitors in neuromorphic circuits [36; 38] is a narrow but highly promising research field. In fact, the possibility of energy-efficient neuromorphic computing with solid-state memcapacitive structures has been demonstrated recently [9].

The mechanical aspect of biological neurons is also becoming increasingly recognized as important in the study of neuroscience. They play a role in several physiological processes, including the generation of action potentials [29]. The mechanical changes in the cell membrane can influence the opening and closing of ion channels, which can impact the electrical signaling of the neuron. The soliton model combines mechanical and electrical factors of signal propagation through axons and has been able to account for some experimental effects in anesthesiology beyond the Hodgkin-Huxley model [2; 3; 16; 26].

In this work, in part inspired by biological neurons, we propose a simple neuromorphic electromechanical model based on a leaky memcapacitor that is capable of achieving the periodic generation of spikes under constant stimulation. The presence of additional --mechanical-- degrees of freedom makes it easier to realize neuromorphic behaviors in artificial electromechanical systems. Applying the methods of dynamical systems, we explore the qualitative features of this model in terms of its fixed points, bifurcations, and limit cycles under DC stimulation. In addition, we demonstrate that the system is capable of complex dynamical adaptations, such as synchronization with periodic external drives, spiking frequency drift, bursting, and other adaptations (for some behaviors, an additional memory circuit element is required).

## II The model

In this section, we introduce a leaky memcapacitor, its model, and the circuit we use to simulate the spiking neurons. The circuit, physical diagram of leaky memcapacitor, and its equivalent circuit are shown in Fig. 1.

The central element of our spiking neuron is a _leaky memcapacitor_ - a capacitor with a plate that moves in response to the force exerted by the internal electric field and restoring force of the spring (the spring constant is \(k\)). In the absence of charges on the plates, \(q=0\), the distance between two plates is \(d\). It is assumed that the displacement of the top plate, \(x\), is positive for the displacement towards the bottom plate, see Fig. 1(a). Thus, the distance between the plates is \(d-x\).

We assume that the capacitor is _leaky_: there is a finite resistance \(R(x)\), which depends on the distance between the plates, dropping rapidly in the 'contact' region, \(x>x_{c}\). The leaky memcapacitor has the _memory_ (which justifies the

Figure 1: (a) An electromechanical leaky memcapacitor connected to a voltage source through a resistor. (b) The equivalent electronic circuit of the leaky memcapacitor: a memcapacitive system, \(\mathrm{C_{M}}\), and memristive system, \(\mathrm{R_{M}}\), connected in parallel. Both systems depend on the same internal state variable \(x\).
term memcapacitor) because the current position of the moving plate (and hence the capacitance itself) depends on the prior history of the charge on the capacitor, which in turn depends on the history of applied voltage. It is this history dependence, the memory, that is responsible for the emergence of complex behaviors that the simple circuit in Fig.1(a) exhibits.

The main operating principle of this device is straightforward. Starting from the equilibrium position at \(x=0\), upon turning on voltage \(V\), the capacitor begins to charge up. This leads to the attraction between its plates, which brings them closer together. For large enough applied voltage, the top plate reaches the contact region, the capacitor discharges, and the plate recoils back from the contact region. Note that the process of charging is limited by the resistor with the resistance \(r\) (Fig. 1(a)).

We found that the proximity-induced attraction between the two surfaces (plates) is helpful to avoid the plates finding a new static equilibrium position in the contact region. We model this part of interaction with a Lennard-Jones-like potential. The short-range attraction between the plates allows for a more effective discharge process, leading to a periodic approach of the top plate to the bottom plate followed by recoil. The exact form of \(R(x)\) and the form of the potential, to a certain extent, do not matter. In the biological context, our memcapacitor may represent two lipid monolayers forming a cell membrane. As the membrane swells or thins, the resistance of the membrane (via ion channels) is affected, as modeled by \(R(x)\).

In our modeling of membrane dynamics, we neglect the kinematic mass (in a biological setting that would be because membranes reside in an aqueous environment that has considerable viscosity). This overdamped regime makes it nontrivial to have oscillatory behavior, and new types of mechanisms are required to perform the periodic spiking that we see possible.

Eqs. (1) and (2) provide the general framework for the description of leaky memcapacitor. For the capacitive and resistive responses, Eq. (1) is written as

\[q =\frac{\epsilon A}{d-x}V_{C}\equiv C(x)\cdot V_{C}\;\;, \tag{3}\] \[V_{C} =\left[R_{m}\cdot\left(\frac{1}{\pi}\mathrm{arctan}\beta\left( x_{c}-x\right)+0.5\right)+\rho_{0}\frac{d-x}{A}\right]I_{M}\equiv R(x)\cdot I _{M}\;\;, \tag{4}\]

where \(q\) is the memcapacitor charge, \(V_{C}\) is the voltage across the leaky memcapacitor, \(I_{M}\) is the leakage current, \(\varepsilon\) is the permittivity, \(A\) is the plate area, and \(R_{m}\), \(\beta\), \(\rho_{0}\) are the parameters defining the memristance. In the normal regime, \(R(x)\approx R_{m}\), while in the contact regime, \(R(x)\) is mainly determined by the \(\rho_{0}\) term; \(\beta\) describes the rate of change of the \(R_{m}\) component in \(R(x)\) in the vicinity of \(x_{c}\).

The leaky memcapacitor is described by a single internal state variable \(x\), which is the displacement of the top plate from its equilibrium position (at \(q=0\)) in the downward direction (see Fig. 1(a)). Its dynamics (corresponding to Eq. (2)) is represented by

\[\gamma\dot{x}=\frac{q^{2}}{\epsilon A}-\frac{\mathrm{d}U(x)}{\mathrm{d}x}\;\;, \tag{5}\]

Figure 2: (a) Resistance as a function of the displacement of the top plate (Eq. (4)). (b) Potential energy as a function of the displacement of the top membrane (Eq. (6)). The insets present a zoomed-in contact region. The dashed line refers to \(x_{c}\).
where \(\gamma\) is the dissipation coefficient, and the potential \(U(x)\) is chosen as

\[U(x)=\frac{1}{2}kx^{2}+4\epsilon_{l}\left[\left(\frac{\sigma}{d-x}\right)^{12}- \left(\frac{\sigma}{d-x}\right)^{6}\right]\ . \tag{6}\]

Here, the first term is the spring potential energy, while the second one is the Lennard-Jones-like potential that we use to describe the contact interaction between the plates. In Eq. (6), \(k\) is the spring constant, \(\epsilon_{l}\) is the depth of the Lennard-Jones potential well, and \(\sigma\) is the distance at which the Lennard-Jones potential energy is zero.

In the following, we measure the distances in the units of \(\sigma\), resistances in the units of \(R_{m}\), \(U\) in the units of \(\epsilon_{l}\), \(k\) in the units of \(\epsilon_{l}/\sigma^{2}\), time in the units of \(R_{m}A\epsilon/\sigma\), charge in the units of \(\sqrt{\epsilon\epsilon_{l}A/\sigma}\), \(\gamma\) in the units of \(\epsilon\epsilon_{l}AR_{m}/\sigma^{3}\), voltage in the units of \(\sigma\epsilon_{l}/\epsilon A\), and current is the units of \(\sqrt{\sigma\epsilon_{l}/\epsilon A}/R_{m}\). For thus defined dimensionless variables and parameters, we use the original notation to minimize clutter in the text.

The parameters used in our simulations are given in Table 1. Fig. 2(a) shows \(R(x)\) as defined by Eq. (4); Eq. (6) is presented in Fig. 2(b).

Note that qualitatively similar results (to the ones presented in this paper) may be obtained using a different choice of parameters and functional dependencies. In particular, we have verified that the results remain nearly the same when the constant resistance for \(x\lesssim 6\) in Fig. 2(a) is replaced with a resistance linearly dependent on \(x\) (for more details, see the Supplemental Information (SI) Appendix A).

To simulate Fig. 1(a) circuit, we use Kirchhoff's voltage law

\[V(t)=r\left(\dot{q}+I_{M}\right)+V_{C}\ \, \tag{7}\]

which is applied together with the equations defining the leaky memcapacitor, Eqs. (3)-(6). The trajectories \((x(t),q(t))\) are found by numerical integration of Eqs. (5) and (7). The integration was performed using the ODE solver ode45 for nonstiff differential equations in MATLAB (version R2022a).

## III Spiking behavior and phase diagram

Figs. 3(a) and (b) show two selected simulation results that demonstrate the transient dynamics from zero initial conditions, \(x_{0}=0\) and \(q_{0}=0\), to the regime of periodic spiking. These results indicate a significant dependence of the spike shape on the magnitude of applied voltage (\(V=8.0829\) in Fig. 3(a) and \(V=15.0111\) in Fig. 3(b)). According to Figs. 3(a) and (b), the spikes are smoother at \(V=8.0829\), and sharper at \(V=15.0111\). The sharper spikes show a close resemblance to biological spikes.

A notable feature of the transient dynamics is the initial sharp increase in the voltage across the leaky memcapacitor, \(V_{C}\). This property (clearly observed in Figs. 3(a) and (b)) is associated with the smaller capacitance at short times (close to \(t=0\)) due to the initially large plate separation. Within the transient region, the leaky memcapacitor adapts to the applied voltage: the Coulomb attraction reduces the distance between the plates what increases the capacitance. The voltage plateau in Fig. 3(a) is close to the bifurcation point at \(V_{l}^{\prime}=7.9582\) (for more information, see SI Appendix A). This explains the relatively long duration of the transient region.

Moreover, in Fig. 3(c) we present an example of non-spiking behavior. In this case, the trajectory ends in a sink (attractor).

We systematically analyzed the behavior of the circuit by studying its fixed points and limit cycles. For this purpose, we used vector field diagrams of solutions and Jacobi matrices (for more information, see SI Appendix A). This analysis has resulted in the phase diagram presented in Fig. 4(a). The diagram indicates the presence of a global limit cycle in a wide range of the applied voltage, from \(V_{1}^{\prime}\) to \(V_{2}\), which accounts for the presence of the spiking regime.

The main features of the phase diagram are as follows. First of all, since the behavior is symmetric with respect to the sign of \(V\), we show only the positive voltage region of the phase diagram. Around \(V=0\), the only global attractor is a sink. As \(V\) increases, at \(V_{0}\approx 3.4156\) there is a bifurcation that nucleates a saddle and spiral source; they do not influence the global attractor, however, until they separate sufficiently at \(V_{1}\) where the saddle splits

\begin{table}
\begin{tabular}{c c c c} \hline Parameter & Value & Parameter & Value \\ \hline \(d\) & 8 & \(r\) & \(10^{-3}\) \\ \(x_{c}\) & 6.4 & \(\rho_{0}\) & \(1.25\times 10^{-4}\) \\ \(\beta\) & \(5\times 10^{4}\) & \(\gamma\) & \(1.25\times 10^{-4}\) \\ \(k\) & 5/6 & & \\ \hline \end{tabular}
\end{table}
Table 1: Parameters used in simulations.
the phase space into two disconnected regions, one of which hosts a limit cycle and the other the original sink that corresponds to a static state, as shown in Fig. 4(b). In SI Appendix B we show the possibility of switching between these two attractors using voltage pulses. Another bifurcation occurs when the sink and saddle point annihilate at \(V_{1}^{\prime}=7.9582\), transforming the limit cycle into the global attractor with stable spikes, as shown in Fig. 4(c). At \(V_{2}\), another bifurcation generating a sink-saddle pair cuts off the limit cycle, shifting the global attractor to the sink, as shown in Fig. 4(d). As \(V\) continuously increases, the saddle and spiral source point move towards each other and then annihilate, leaving the sink alone (for more information, see SI Appendix A).

To better understand the properties of spikes, the _natural_ frequency of spikes, \(w_{natural}\), was calculated as a function of the applied voltage (see Fig. 5(a)). Interestingly, the calculated points are distributed in the half-of-the-oval shape in the frequency-voltage plot. Fig. 5(a) shows that the frequency approaches zero when \(V\to V_{1},V_{2}\). The Fourier transform of the voltage across the memcapacitor is presented in Fig. 5(b). Qualitatively, the whole spiking regime can be divided into three parts, I, II, and III, that are different by the pattern of spikes (see Figs. 5(c)-(e)). The "negative spikes" regime, I, and the "positive spikes" regime, III, are connected by the regime of more symmetric (harmonic) spikes, II.

## IV Synchronization with external source

To study the synchronization with an external source, an ac voltage was added to the constant driving voltage \(V_{dc}\), \(V(t)=V_{dc}+\delta V\sin(\omega_{source}t)\), with \(\delta V=0.1155\). To initialize the system close to the limit cycle, we used the initial conditions \(x_{0}=6.60\) and \(q_{0}=2.02\). The Fourier transforms for the regimes I-III of oscillations in Fig. 5(a) are presented in Fig. 6. When the circuit is in regime II (Fig. 6(b)), i.e. away from the thresholds, the synchronization occurs only when the source frequency is very close to the spike frequency, \(\omega_{source}\approx\omega_{natural}(V_{dc})\). In this regime, the Fourier plot features various combinations of integer multiples of \(\omega_{source}\) and \(\omega_{natural}\), \(N\omega_{natural}\pm M\omega_{source}\), where \(M,N=0,1,2,...\)

In the other two cases in Fig. 6, the external source has a much stronger influence on the spike generation. Figs. 6(a) and (c) represent the cases when \(V_{dc}\) is close to the thresholds \(V_{1}\) and \(V_{2}\). In these cases, the regions of synchronized driving frequencies are largely extended in both integer multiples and rational fractions of the self-oscillation frequencies. In regime I, when the driving frequency is small, the system responds harmonically to the drive, as shown in Fig. 6(a), where only the component of \(\omega_{source}\) is evident. This is because once \(V\) drops below \(V_{1}\), the spiral source and the saddle vanish simultaneously, and the system can easily switch to the sink, which is away from the contact regime as depicted in Figs. 4(a) and (b), and may not return. In contrast, in regime III, when \(V\) crosses over \(V_{2}\), the bifurcation of a sink-saddle pair appears in the contact regime on the previous limit cycle, as depicted in Figs. 4(a) and (d), thus the system will not go away from the contact regime, which results in spiking waveform even at low source frequency. Further detailed Fourier transforms and relevant analyses of these three regimes are discussed in SI Appendix C.

Figure 3: The response to a step-like voltage applied at \(t=0\): starting from zero initial conditions (\(x=0\), \(q=0\)), the circuit transits to (a), (b) a periodic spiking regime or (c) static regime. These plots were obtained using (a) \(V=8.0829\), (b) \(V=15.0111\), and (c) \(V=7.8520\). In the top pannels, The black dashed line refers to \(d\) and the brown one refers to \(x_{c}\).
## V Other types of dynamics

In this section, we show that the dynamics of Fig. 1(a) circuit can be further enriched by the use of additional components with memory. Indeed, the replacement of the resistor \(r\) in Fig. 1(a) with a memristor may completely change the pattern of spikes to bursting (a pattern of firing wherein the periods of rapid spiking are separated by quiescent periods). Below, we introduce two memristive behaviors (without and with threshold with respect to \(I_{r}\)) and study their effect on the spike generation.

In principle, it is evident that at a suitable constant applied voltage \(V\), the resistance \(r\) controls the spike generation. As in the limiting cases of \(r\to 0\) and \(r\rightarrow\infty\) the dynamics of Fig. 1(a) circuit should not be oscillatory, one can assume that the spiking behavior occurs within certain resistance thresholds, say, when \(r_{2}<r<r_{1}\). Consequently, a suitable memristor (e.g., with resistance changing across \(r_{1}\) and/or \(r_{2}\)) may be used to control the pattern of spikes.

For the sake of simplicity, in what follows the resistance is used as the internal state variable of memristor (such models are known in the literature [37]). In the first model, the memristor dynamics is defined by the equation

\[\dot{r}(t)=-\alpha_{1}I_{r}(t)+\lambda_{1}\int\limits_{0}^{t}e^{-\gamma(t- \tau)}\left(r_{0}-r(\tau)\right)\mathrm{d}\tau\;\;, \tag{8}\]

where \(\alpha_{1}\) is the current proportionality coefficient (can be of any sign), \(\lambda_{1}\) is the (positive) relaxation coefficient, \(\exp\left(-\gamma\left(t-\tau\right)\right)\) is the memory kernel, \(\gamma\) is non-negative constant, \(r_{0}\) is the equilibrium resistance, and \(I_{r}(t)=\dot{q}+I_{M}\) is the current flowing through \(r\). It is assumed that \(r\) is confined to the interval from \(r^{\prime}\) to infinity, where \(r^{\prime}\) is the

Figure 4: (a) An overall picture: the fixed points and attractors depending on the applied voltage \(V\). The phase diagrams at (b) \(V=7.0437\) (the sink and limit cycle case in (a)), (c) \(V=7.9674\) (the limit cycle case in (a)) and (d) \(V=15.5000\) (the sink case on the right of the limit cycle one in (a)). In (b) the phase space is divided into two parts by the flow lines towards the saddle, as depicted with a dashed line. The semi-transparent gray thick line represents the limit cycle. The saddle, the sink, and the spiral source are labeled by black arrows, and the yellow star inside the limit cycle is located at the spiral source. The blue arrows in (b) and the orange arrows in (d) depict the shift direction of the three fixed points as \(V\) increases. In (d), vector fields are added in small blue arrows to help show the fixed points. The initial conditions of the solutions are set discretely along the edge, as well as around the spiral source. In (b)-(d), the evolution time \(t\) was \(0.05\).
minimal value of \(r\). In the right-hand side of Eq. (8), the first term causes the drift of \(r\) mostly due to the spikes. The second term incorporates the relaxation to \(r_{0}\) and memory on prior states.

Figs. 7(a) and (b) show results of simulations with Eq. (8) model for \(r\). Roughly the picture is as follows: The spikes starting at \(t=0\) cause \(r\) to decrease. As soon as \(r<r_{2}\), the spiking terminates and \(r\) starts relaxing towards \(r_{0}\) with a delay determined by \(\gamma\). When \(r\) crosses \(r_{2}\) again (in the opposite direction), the circuit transits back into the spiking mode. This process repeats periodically.

Figure 5: (a) Spike frequency as a function of \(V\), and (b) Fourier transform of \(V_{C}\) as a function of \(V\). The dashed lines, from left to right, refer to \(V_{1}\), \(V_{1}^{\prime}\) and \(V_{2}\), respectively. Here, \(P_{1}\) is the single-sided amplitude of the Fourier transform. In these calculations, to keep the system close to the limit cycle, the initial condition was selected as \(x_{0}=6.6000\) and \(q_{0}=2.0207\). The evolution time \(t\) was \(1.5\), and we skipped the initial transient interval. (c)-(e) Steady-state oscillations at (c) \(V=7.0183\) (regime I), (d) \(V=11.5470\) (regime II) and (e) \(V=15.0561\) (regime III). In (c)-(e), the black dashed line refers to \(d\) and the brown one refers to \(x_{c}\).

Figure 6: Synchronization with external source: The Fourier transform of \(V_{C}\) at (a) \(V_{dc}=7.0645\) (regime I), (b) \(V_{dc}=11.5470\) (regime II) and (c) \(V_{dc}=15.0688\) (regime III). \(P_{1}\) is the single-sided amplitude of Fourier transform. In these calculations, to keep the system close to the limit cycle, the initial condition was selected as \(x_{0}=6.6000\) and \(q_{0}=2.0207\). These plots were obtained using the integration time \(t\) of \(2\); the initial interval of transient dynamics was omitted. In (a)-(c), the horizontal and vertical dashed lines correspond to \(\omega_{natural}(V_{dc})\).
Next, we consider the second model for \(r(t)\). A deeper analysis of Fig. 1(a) circuit (with \(r\) as a control parameter) has revealed an additional threshold \(r_{1}^{\prime}\) between \(r_{2}\) and \(r_{1}\) such that there are two attractors (a sink and limit cycle) for \(r_{1}>r>r_{1}^{\prime}\). In the limit cycle, certain variables such as \(I_{r}(t)\) have higher values for the duration of one spiking period compared to those in the sink. Therefore, we can set a threshold, \(I^{\prime}\), to induce the cycling between them. If

Figure 7: Achieving complex spiking behaviors by introducing memory into the resistor. (a), (b) Method of Eq. (8). (c), (d) Method of Eq. (9). Panels (b) and (d) are zoom-ins of (a) and (c), respectively. In the inset of (d), the dashed line refers to \(I^{\prime}\), and during the short plateau, the limit cycle has just disappeared and the system is switching to the sink. In the simulation (a), the parameters were set as \(V=13.8564\), \(\alpha_{1}=3.4641\times 10^{-6}\), \(\lambda_{1}=1.600\times 10^{5}\), \(r_{0}=10^{-3}\), \(\gamma=0\) and \(r^{\prime}=0.8r_{0}\). In the simulation (c), \((\arctan\lambda^{\prime}(I_{r}(t)-I^{\prime}))/\pi+0.5\) is used instead of step functions in Eq. (9), and the parameters were set as \(V=7.9674\), \(I^{\prime}=10.3923\), \(\lambda^{\prime}=1000\), \(\alpha_{2}=1.6000\times 10^{4}\) and \(\lambda_{2}=2000\). Initial conditions are both set \(x=0\), \(q=0\) and \(r_{0}=10^{-3}\).
\(I_{r}(t)>I^{\prime}\), we can set \(\dot{r}>0\), otherwise \(\dot{r}<0\), like

\[\dot{r}(t)=\begin{cases}\alpha_{2}r(t)^{2},&I_{r}(t)>I^{\prime}\\ -\lambda_{2}r(t),&I_{r}(t)<I^{\prime}\\ 0,&I_{r}(t)=I^{\prime}\end{cases}\, \tag{9}\]

where \(\alpha_{2},\lambda_{2}\) are positive coefficients. It is assumed that \(r\) is confined to the interval from \(r^{\prime}\) to infinity, where \(r^{\prime}\) is the minimal value of \(r\). This mechanism allows \(r(t)\) to increase from \(\{r^{\prime}_{1},r_{1}\}\) to above \(r_{1}\), switching the system from the limit cycle to the sink; after a while \(r(t)\) begins to decrease, closing the cycle. The simulations of this scheme for \(\dot{r}\) are shown in Fig. 7(c) and (d). Note that the second method works only when \(R_{m}\) is finite (for more details, see SI Appendix A).

## VI Conclusion

In summary, we have proposed a leaky memcapacitor - an electromechanical crossbreed of memcapacitor and memristor - that can generate neuromorphic spikes. Its model is based on the potential that combines linear elasticity with nonlinear Lennard-Jones-like interaction between the plates at short distances, attempting to represent realistic interaction potential. Thanks to the presence of the non-linear interaction, the dynamical behavior of the system in the contact region is different from the one when the plates are relatively far from each other. This helps achieve a stable spiking behavior when a constant voltage is applied.

In order to thoroughly understand the spiking behavior, we have conducted the stability analysis in the \((x,q)\)-space and discovered several interesting regimes characterized by different configurations of fixed points and attractors. We have shown that for some ranges of parameters, one can use a voltage pulse to switch the system from a sink to a limit cycle and vice-versa. We have also found that the spike shape depends on the applied voltage.

An important feature of the system is that the spike frequency may adapt to the external perturbation frequency (depending on the model and excitation parameters). A rich dynamical behavior has been observed including synchronization when a small amplitude ac signal was added to the constant driving voltage. In addition, replacing constant external resistor by a memristor extends the variety of spike waveforms that the circuit generates. With this modification, the circuit can be tuned to mimic behaviors of some types of biological neurons.

Overall, the system introduced in this article provides a new avenue for the practical realization of neuromorphic devices based on memcapacitive and memristive effects. Our study may lead to novel energy-efficient realizations of neural dynamics with electromechanical structures, including artificial analogs of biological membranes.

## Conflict of Interest Statement

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.

## Author Contributions

All authors contributed to conceiving the idea, formulating the model, its qualitative analysis, and writing the manuscript. Z. Z. has performed all numerical simulations.

## Funding

I. M. acknowledges funding from the Materials Sciences and Engineering Division, Basic Energy Sciences, Office of Science, US DOE.
## Data Availability Statement

The Supplemental Information contains additional data supporting this manuscript. Further inquiries can be directed to the corresponding authors.
