# Dynamics of neural motifs realized with a minimal memristive neurosynaptic unit

Jiaming Wu, Adrien D'hollande, Haoran Du, Marcelo Rozenberg

► To cite this version:

Jiaming Wu, Adrien D'hollande, Haoran Du, Marcelo Rozenberg. Dynamics of neural motifs realized with a minimal memristive neurosynaptic unit. Physical Review Applied, 2025, 23 (3), pp.034030. 10.1103/PhysRevApplied.23.034030. [hal-05345008]

HAL Id: hal-05345008

https://hal.science/hal-05345008v1

Submitted on 3 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Dynamics of neural motifs realized with a minimal memristive neurosynaptic unit

Jiaming Wu, Adrien d'Hollande, Haoran Du, and Marcelo Rozenberg

Université Paris-Saclay, CNRS, Laboratoire de Physique des Solides, Orsay 91405, France

<figure><img src="image_2.png" /></figure>

(Received 5 July 2024; revised 9 January 2025; accepted 13 February 2025; published 14 March 2025)

This paper is a contribution to the Physical Review Applied collection titled Physics-Inspired Computing.

Spiking neural networks (SNNs) require both neurons and synapses. Here we introduce a synaptic unit with both excitatory and inhibitory biomimetic dynamics (Rall's function). Combined with our previously introduced memristive spiking neuron, they complete a minimal model platform to build general SNNs. Our hardware neurosynaptic (NS) unit is fully analogue and implements a textbook theoretical model. It avoids the mismatch problem and is easily tunable at biocompatible timescales. We demonstrate the modularity of our NS-unit platform by building various basic neurocomputing primitives produced by basic motifs and their combination. With our novel hybrid theoretical and experimental approach we make three notable contributions: (1) We demonstrate that recursive inhibition (adaptation) and excitation can both be formulated as self-consistent dynamical problems. (2) We demonstrate that a minimal dynamical memory can be implemented with just one recursive single neuron. (3) We show that adopting a biomimetic membrane current in adaptive neuron models reveals an unknown connection to a cornerstone of physics, the harmonic oscillator model. Our NS unit realizes a building block for midsize SNNs of arbitrary geometry. Its simple design, along with the wide accessibility of ordinary electronic components, make our methodology an attractive platform to build neural interfaces for biological neuroscience, medical devices including deep brain stimulation, robotics, and artificial intelligence systems such as reservoir computing. We provide a bill of materials and printed-circuit-board designs to implement the device.

DOI: 10.1103/PhysRevApplied.23.034030

## I. INTRODUCTION: NEUROCOMPUTING PARADIGMS AND THE PROBLEM OF COMPUTING WITH SPIKES

Neuromorphic computing (NC) is a field aimed at implementing systems that can perform functions taking inspiration from the nervous system and brain. Examples range from face recognition through autonomous navigation to robotic motion, and many more that are revolutionizing modern life. To implement those systems, there are two main approaches. One is based on software running on digital computers, and the other is fabricating neuromorphic hardware that mimics the components of a neural network. In both cases, NC systems rely on models from theoretical neuroscience, computer science, or physics. Every model has a definition of its components, called neurons and synapses. Neurons are the degrees of freedom of a neural network, since they can be in different states. Synapses describe the couplings and interactions between them (see Fig. 1).

We recall a paradigmatic example of an NC system, the Hopfield model, to clarify the roles of neurons and synapses, which are often a matter of confusion. This is a model for associative memory, introduced in the 1980s

[1]. It bears close similarity with magnetic spin models in physics (see Fig. 1). The neurons are defined as discrete variables, where $S_i = \pm 1$ denotes the $i$th member of the network. The synapses are the interactions between neurons, where the real parameter $J_{ij}$ denotes the interaction between neurons $i$ and $j$. The model has an energy given by the expression $E = (1/N) \sum J_{ij} S_i S_j$, where $N$ is the number of neurons. Each memory corresponds to a local minimum of the energy landscape. The genius of Hopfield was to formulate a simple recipe to choose the $J_{ij}$ couplings that provide such a landscape. The model can recall a memory by association. Namely, given an arbitrary initial state $S_i$ for $i = 1, ..., N$, an algorithm of energy minimization will retrieve its associated memory.

### A. Neurocomputing with software

The Hopfield model, as any of the modern deep neural networks for artificial intelligence, is an algorithm implemented in software that runs on digital computers. The computers may be conventional or the so-called “neuromorphic chips,” such as TrueNorth or Loihi, which are digital processors with an optimized architecture to run those types of models [2]. Bitlike variables and $J_{ij}$ couplings are the simplest modeling for neurons and synapses and allow for computation with massive neural networks, such as those used in artificial intelligence. However,

<footnote>*Contact author: Marcelo Rozenberg, mjrozenberg@gmail.com, marcelo.rozenberg@universite-paris-saclay.fr</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. Schematic view of the presynaptic (red) and postsynaptic (blue) neurons. The neurons are coupled by a synapse. The action potential reaching the end of the axon terminals of the preneuron induces a synaptic input current in the postneuron, which can be excitatory (positive) or inhibitory (negative). The inset shows the analogy with the Hopfield model, where neurons are spins $S_i$, and synapses are the couplings $J_{ij}$.</figcaption></figure>

severe limitations emerge quickly if one wants to adopt a more biomimetic spiking-neural-network (SNN) model.

In fact, neurons and synapses are qualitatively different than the simplified description of the Hopfield model: they are dynamical systems. Neurons are a special type of cell in the body that are characterized by having electrical activity. They are not just active or inactive, but emit action potentials that are fast voltage spikes. Synapses involve the dynamical process of transferring neurotransmitters between a presynaptic and a postsynaptic neuron, so they are not just real-valued parameters (see Fig. 1) [3].

The dynamical models of spiking neurons range from schematic to biorealistic. In the first case, the most paradigmatic one is the leaky-integrate-and-fire (LIF) model [3,4]. For the latter, the most famous one is the Hodgkin-Huxley model, which describes the generation of action potential in the giant axon of the squid [5]. These models are described by multiple nonlinear differential equations. Similarly, to describe the synaptic couplings, further differential equations are needed. Hence, to numerically study an SNN requires the solution of a challenging system of multiple coupled nonlinear differential equations. The scaling of these SNNs is typically bad, as the number of differential equations grows rapidly with the size of the network and its synaptic connectivity. Possibly an even more serious problem is the scaling of the simulation time. An action-potential spike typically lasts 1 ms, and the firing rate of a neuron is roughly between 1 and 100 Hz. Thus a time step of 0.1 ms is often adopted. The timescale associated with the global behavior can be considered of the order of minutes, so there is at least a 10⁶ factor for the number of simulation time steps, just for one neuron. However, this scaling quickly worsens with a growing number of neurons, at least by a factor of N, in the most favorable case. Nevertheless, there are several implementations

of useful numerical simulators of SNNs available (see [6] for a recent review and benchmark).

It is not easy to find out what is the largest SNN size that current simulators can handle [6], as this depends on the model, the connectivity, the computer architecture, the power, etc. Whatever this estimate may be, it would become severely reduced if one considers continuous models, without hard resets of variables, where the spike emission is not described by time stamps and synapses are dynamical. Namely, where the full nonlinearity of the problem is properly treated.

Finally, another severe limitation imposed on software approaches is the numerical and electric power requirements needed to carry out the computation. Power efficiency becomes a relevant aspect, especially for practical applications such as navigation control, neuroprostheses, wearable implants, brain-machine interfaces, etc. [7,8].

## B. Neurocomputing in hardware

A very different approach to NC is based on the idea of implementing neurons and synapses directly in hardware. In fact, a hardware implementation has the very attractive feature of *ideal scaling* because time represents itself. In other words, the system evolves in real time and the time evolution occurs irrespective of the number of neurons in the network. If it takes one minute to simulate the behavior of one neuron, it takes the same one minute to simulate the behavior of a network of size N.

Modern silicon electronic technology based on complementary metal oxide semiconductor (CMOS) systems can implement, *a priori*, almost any desired mathematical model. This approach has had several decades of development, starting with the pioneering work of Mead, who coined the term “neuromorphic electronics” [9]. He observed that transistors operating in the subthreshold regime, i.e., not as digital switches but as variable conductances, show a behavior analogous to the activation of ionic channels in neuron membrane models. There are currently many electronic chips that integrate millions of *silicon neurons*, which typically require tens of transistors each [10,11]. Surprisingly, the implementation of dynamical synapses in hardware has received comparatively less attention [12], although they are a key feature in theoretical neuroscience [13].

A remarkable example of this type of hardware implementation is the DYNAP-SE chip, which counts 4 × 256 neurosynaptic computing units [14]. Each one implements an adaptive-exponential integrate-and-fire neuron [3] that generates discrete spike events, plus a dynamical synapse that can implement biomimetic currents. This chip exhibits the appealing feature of ultralow power consumption. In contrast, some significant inconveniences of this approach are that it requires a high level of technical expertise and of fabrication capabilities, beyond the reach of most
research labs, and that it lacks flexibility, as it takes years to design and implement each electronic chip. Nevertheless, the most significant of all issues, common to all other electronic implementations that work in the transistor subthreshold regime, is that it suffers from a severe problem of device mismatch [15].

Here we propose to take a drastically alternative approach to fabricating, studying, and understanding SNNs. While current CMOS VLSI may be considered as a *top-down* method, where fabricating a large number of units is a main goal, here, we propose a *bottom-up* approach. Our methodology is constructive, with neurosynaptic units that can be modularly interconnected to implement neurocomputing primitives. Crucially, we can do this in a fully controlled manner, including the homogeneity between the spiking units, hence mitigating the mismatch problem. Thus, our approach to implementing SNNs in hardware qualitatively differs from all others in multiple aspects: it is of extreme simplicity, extreme flexibility, extreme tunability, controlled reproducibility, very low fabrication cost, and based on theoretical neuroscience textbook modeling.

Here, we shall introduce and demonstrate our hybrid theoretical and experimental methodology by implementing in hardware all the basic neurocomputing motifs from theoretical neuroscience. The originality and value of our approach will be highlighted by the relevant novel results and insights that emerge from the study of dynamical systems.

Specifically, here we shall demonstrate the following important results.

(1) We demonstrate that the minimal recursive SNN, i.e., a single neuron with either autoinhibition (adaptation) or autoexcitation, can be mapped onto a self-consistent dynamical problem with a geometrical representation. Based on this formulation and its hardware implementation we demonstrate the two following important results.

(2) For the case of autoexcitation, we obtained the unprecedented realization of a single-neuron stable dynamical memory (also known as a working memory). This is significant, because in the current literature it is widely assumed that, to achieve a stable working memory, a recursive neuron assembly is needed [3]. This finding may also explain the recent striking experimental observation of the bump of activity in the head direction system of *Drosophila*. This is realized in a neural structure called the ellipsoidal body, which counts with only a handful of neurons [16].

(3) For the case of autoinhibition, we made two significant findings. First, we show that our model is a hardware counterpart of paradigmatic models of theoretical neuroscience, the Izhikevich and the AdEx [3,17]. Similarly to those models, the adaptation is mediated by a spike-dependent membrane current. But, improving on them, we model the current with the more biorealistic Rall's function

[13]. This yields the second significant finding, where we theoretically show that the model can be mapped to a cornerstone model of physics, namely, the driven harmonic oscillator with friction. We also demonstrate this experimentally, by the measured currents in our hardware circuits that realize an intrinsic burster or pacemaker state.

(4) We show how the single-neuron motifs can be leveraged as components of higher-level two-neuron motifs, which produce three different types of central pattern generators (CPG). Our results demonstrate two features: (i) the modular constructive nature of our methodology; and (ii) that our CPGs can easily work at biocompatible timescales, making them ready for deployment in practical applications.

The paper is organized as follows. In Sec. II we first describe the neurosynaptic unit model and its circuit implementation, emphasizing the connection to mathematical models and its full parameter tunability. Then, in Secs. III and IV we demonstrate the implementation of the simplest motifs of spiking neural networks, that realize *neuro-morphic primitives*, i.e., basic neurocomputing functions. Then, in Sec. V we discuss how our present work is related to the intense activity in the search for memristive materials. We also make some general remarks about the positioning of our approach with respect to other hardware in regard to power efficiency. In the final Sec. VI we conclude with a recap on the main contributions of our present work both for science and in neuromorphic electronics. Finally, in the Appendix we provide technical details on all our circuits and a list of materials, which are also available upon reasonable request.

## II. NEUROSYNAPTIC UNIT: THEORETICAL MODEL AND ITS PHYSICAL EMBODIMENT

### A. Theoretical model

The neurosynaptic (NS) hardware unit takes inspiration from a basic model in Wilson's textbook on theoretical neuroscience (cf. Ch. 12 of Ref. [13]). It is a set of differential equations for a spiking neuron of type I supplemented with synaptic currents. A type-I neuron, such as a cortical neuron [13,18], characteristically fires action potential spikes at arbitrarily low rates.

When a presynaptic (or upstream) neuron emits a spike, which is a very fast event, typically lasting 1 ms, the action potential liberates a certain amount of neurotransmitters (see Fig. 1) into the synaptic cleft, which is about 20 nm wide [19]. Those molecules reach the dendrites of the postsynaptic (or downstream) neuron where they open ionic channels that induce synaptic currents integrated in the cell body. The neurotransmitters can be either excitatory or inhibitory, depending on whether they contribute or prevent the emission of a spike by the postsynaptic neuron. This depends on the sign of the synaptic current: if positive, or depolarizing, it brings the potential of
the neuron closer to the firing threshold; if negative, or hyperpolarizing, it has the opposite effect.

This synaptic current process typically occurs in the range of tens to hundreds of milliseconds, thus at a much longer timescale than the spike emission, and has its own dynamics. This introduces an additional and relevant timescale in SNN dynamics, which is well recognized to have significant consequences [20].

In the 1960s Rall introduced a set of differential equations to model this phenomenon, which has been widely adopted [3,13,21,22]. When applied to a delta function representing the sharp spike emission, they provide as a solution a synaptic current with the functional form known as Rall's alpha function, $\alpha(t) = (t/\tau_s^2) \exp(-t/\tau_s)$. Thus, the sharp (~ 1 ms) spike emitted by the preneuron is perceived by the postneuron as a broad synaptic current input that peaks and decays in a typical timescale $\tau_s$ (~ 10–100 ms).

Rall's differential equations for the synaptic current are incorporated in Wilson's book model (cf. Eq. (12.18) of Ref. [13]). Those equations, along with the expression for a general integrate-and-fire spiking neuron from Gerstner's book (cf. Eq. (5.2) of Ref. [3]), conform to the mathematical model that underlies our hardware implementation. The set of equations reads

$$C_m \frac{dV_m}{dt} = f(V_m) + I_\alpha + I_0, \quad (1)$$

$$\tau_s \frac{dI_e}{dt} = -I_e + I_\delta H(V_{\text{pre}} - \Omega), \quad (2)$$

$$\tau_s \frac{dI_\alpha}{dt} = -I_\alpha + I_e. \quad (3)$$

Here $V_m$ is the neuron's membrane potential; $C_m$ is the membrane capacitance; $\tau_s = R_s C_s$ is the synapse time constant, with $C_s$ and $R_s$ being the capacitance and leak resistance of the synapse, respectively; $I_0$ is an external input current to the neuron; and $I_e$ and $I_\alpha$ are dynamical synaptic currents (analogous to $f$ and $g$ in Wilson's book [13]). Also closely following the notation of Wilson, $H$ denotes a step function (or Heaviside function) that is unity if the presynaptic voltage $V_{\text{pre}}$ is larger than the synaptic threshold $\Omega$, or zero otherwise [13]. Thus, the presynaptic action potential is modeled as a narrow pulse of rectangular shape that represents a $\delta(t)$-like spike. The condition ($V_{\text{pre}} > \Omega$) sets the width of the narrow pulse, which we call $\tau_a$, and its magnitude we denote by $I_\delta$. Hence, this spike carries a charge

$$Q_N = I_\delta \tau_a \quad (4)$$

that physically may represent the amount of neurotransmitter liberated by the preneuron into the synapse.

The current $I_e$ results from the first leaky integration equation (2) of the sharp rectangular current pulse. Hence,

with $\tau_a \ll \tau_s$, one obtains the exponential current

$$I_e(t) = \frac{Q_N}{\tau_s} \exp(-t/\tau_s), \quad (5)$$

This current is the input to the third equation (3), so after this second leaky integration one gets the alpha function [3,13,22]

$$I_\alpha(t) = \frac{Q_N}{\tau_s} \frac{t}{\tau_s} \exp(-t/\tau_s). \quad (6)$$

These currents have biomimetic line shapes [3,13]. They realize the synaptic currents of four basic types, namely, fast exponential excitatory (AMPA) or inhibitory (GABA$_a$), and slow alpha-function excitatory (NMDA) or inhibitory (GABA$_b$) (cf. Ch. 3 of Ref. [3]).

## B. Analog leaky-integrate-and-fire model

To fully specify the model, we need to define the nonlinear function $f(V_m)$ in Eq. (1) that generates the spikes. In the interest of simplicity, we start by recalling the leaky-integrate-and-fire (LIF) model, which is the most widely used. In the LIF model, shown in Fig. 2, the body (soma) of the neuron is represented as a capacitor $C_m$ that charges by integrating the input ionic currents arriving at its dendrites. This integration suffers leakage losses, which are represented by the membrane resistor $R_m$. The third component of the LIF model is a voltage-dependent switch $S$ that is normally open, but closes when the potential $V_m(t)$ on the capacitor reaches a threshold $V_{\text{th}}$. At that point, one says that a spike is “fired” at a time $t = t_{\text{fire}}$, the voltage $V_m$ is reset to a low resting value, and the switch is set back to open. It is important to realize that the LIF model is a mathematical model and the spike is an abstract “event” defined by the condition $V_m = V_{\text{th}}$ [3]. In other words, there is no dynamical description of the spike [23]. In this model, the nonlinear function $f(V_m)$ is implemented by the V-dependent switch.

Here, we shall adopt the above idealized description as a basis for our hardware implementation, which we shall also augment to provide an embodiment of the physical emission of an action-potential spike (see Fig. 2). The first step is, as shown in Fig. 2, to add a small “load” resistor with $R_a \ll R_m$. Then, when $S$ closes at $V_m = V_{\text{th}}$, the $C_m$ quickly discharges on $R_a$ and the current on this resistor produces a voltage spike. The fast discharge time is $\tau_a = R_a C_m$, which provides a definite parametrization to the width of the $\delta$-like spike that we mentioned above Eq. (4). Interestingly, we may give this spike generation mechanism a biological interpretation as the action potential that is initiated in the axon hillock [24]. The axon hillock is characterized by a high concentration of Na channels, over 100 times more than in the cell body (soma) [25]. We can thus think of the closing of the V-dependent
<figure><img src="image_4.png" /><figcaption>FIG. 2. The analog-LIF model: (a) The leaky-integrate-and-fire model with a voltage-dependent switch S is augmented with a discharge resistance $R_a$ to produce an explicit action potential (left panel). The switch can be implemented by a thyristor T with a resistor $R_m$ between its gate and anode. This resistor also provides the leak resistance of the model, while the thyristor is not conducting (middle panel). The combination of $T + R_m$ realizes the concept of a memristor, with the functionality of a voltage-gated conductance, central to all neuron models (right panel). (b) LIF behavior: The input current pulses (bottom panel) are leaky integrated, producing the increase of the membrane voltage $V_m$ (center panel). When $V_m$ reaches the dynamical threshold $V_{th}$ (vertical arrows), the forward resistance of the T collapses and there is the emission of an output of a continuous action-potential spike $V_{ap}$ (top panel). Circuit technical details and parameter values for all figures are provided in the Appendix.</figcaption></figure>

switch as representing the massive opening of those channels and the enhanced membrane conductance at the axon hillock described by the small resistance $R_a$. Despite the extended use of the LIF model, we are not aware of this simple extension.

The second step to specify the embodiment of the model is to implement the V-dependent switch. The simplest circuit to do this exploits the concept of a volatile memristor, which is a V-dependent resistance (or V-gated conductance) [26]. We have recently, demonstrated that a combination of a thyristor and a resistor can implement such a memristive two-terminal device that provides the desired functionality [27,28].

A voltage-gated conductance is in fact a key ingredient common to all conductance-based neuron models, such as Hodgkin-Huxley, Morris-Lecar, Hindmarsh-Rose, etc. [13]. Those theoretical models are formulated with

differential equations that describe the spike emission as a continuous dynamical process. It results from the strong nonlinearity of the equations, namely the property of excitability [13]. We note that simplified models, such as Izhikevich and AdEx [3], which are popular for numerical simulations, do not fall in that category. This is because they require a discontinuous reset, for the sake of their numerical simplicity [3,23].

Our memristor spiking model implementation is shown in Fig. 2 and it is described by the following equations:

$$C_m \frac{dV_m}{dt} = - \frac{V_m}{R_m[S] + R_a} + I_0, \quad (7)$$

$$V_{ap} = \frac{V_m}{R_m[S] + R_a} R_a. \quad (8)$$

Here $R_m[S]$, with high (hi) or low (lo) S, denotes a volatile memristor whose resistance can switch between two states: a high resistance $R_m[hi]$ and low resistance $R_m[lo]$. The commutation depends on the applied voltage, thus $S = S(V_m)$, implementing a V-dependent conductance [28]. Equation (8) is a voltage divider providing the action potential $V_{ap}$ on the axon-hillock resistor $R_a$. When the memristor is in the $R_m[hi] \gg R_a$ state (i.e., open switch), then $V_{ap}$ is negligible. When it commutes to $R_m[lo] \ll R_a$ (i.e., closed switch), there is an action-potential spike explicitly described by $V_{ap}(t)$ (see Fig. 2).

Notice that the first differential equation above may seem linear, but it is not. In fact, the nonlinearity, necessary for the excitability of the model, is provided by the commutation property of the memristor device [28].

The circuit implementation of the LIF model that we just described features a fully continuous spike generation mechanism. Therefore, our circuit model, although not bio-realistic but biomimetic, can be considered as a member of the class of conductance-based neuron models and we call it the analog-LIF (aLIF) model. In Fig. 2 we show the measured traces of the neuron circuit, which shows the leaky integration of input current pulses and the resulting output emitted spikes $V_{ap}(t)$.

### C. Physical embodiment of the neurosynaptic model equations

One of the main goals of the present work is to provide a simple and transparent physical instantiation of the set of equations (1)–(3), which are the mathematical reference. We would like that each circuit component could be directly related to the theoretical model parameters. Furthermore, we would like the hardware implementation to be modular, such that the subcircuit units can be associated to build arbitrarily spiking neural networks, and eventually a neurocomputer. We shall see in the rest of the paper that all these goals are accomplished by the NS unit that we describe below.
For the spiking generation mechanism, which is encoded in Eq. (7), we adopt the aLIF circuit model discussed in the previous section. Thus, we rewrite below the reference model equations (1)–(3) including the aLIF and renaming some variables for convenience:

$$C_m \frac{dV_m}{dt} = -\frac{V_m}{R_m[S] + R_a} + I_{s2} + I_0, \quad (9)$$

$$V_{\text{post}} = \frac{V_m}{R_m[S] + R_a} R_a, \quad (10)$$

$$\tau_s \frac{dI_{s1}}{dt} = -I_{s1} + I_{\delta}H(V_{\text{pre}} - \Omega), \quad (11)$$

$$\tau_s \frac{dI_{s2}}{dt} = -I_{s2} + I_{s1}. \quad (12)$$

Here $V_{\text{post}} = V_{\text{ap}}$ denotes the output postsynaptic action potential, and $I_{s1}$ and $I_{s2}$ denote the synaptic currents of the two leaky-integration stages.

Similar to what was provided for spike generation via the aLIF circuit, the synaptic current equations (11) and (12) also have a concrete physical embodiment, as we show in the circuit in Fig. 3. We indicate in the figure the correspondence between the equations of the model and the circuit blocks. We shall now describe the implementation in detail. For the sake of clarity, we shall explain the circuit in Fig. 3 from input to output, i.e., from the presynaptic neuron on the left ($V_{\text{pre}}$) to the postsynaptic neuron on the right ($V_{\text{post}}$).

<figure><img src="image_5.png" /><figcaption>FIG. 3. (a) Schematic diagram of the neurosynaptic unit circuit blocks and their relation with the biological neuron. (b) Physical embodiment of model equations (9)–(12) with subcircuit blocks denoting the pre- and postsynaptic neurons and the synaptic block in between. The synaptic block has two stages that perform a leaky integration each. The input presynaptic potential produces a $\delta(t)$-like current spike. A double leaky integration of the pulse provides an $\alpha(t)$-like synaptic current into the postsynaptic neuron, which also receives the external input $I_0$.</figcaption></figure>

## 1. Neurosynaptic circuit blocks

a. The first block. A is an aLIF representing the presynaptic neuron. It generates the action potential spike $V_{\text{pre}}$ that is the input to the NS unit. It was described in detail before.

b. The second block. B implements the differential equation (11) of our model. The current term $I_{s1}(t)$ is the dynamical variable and $I_{\delta}H(V_{\text{pre}} - \Omega)$ is an independent term.

The $I_{\delta}(t)$ currents are rectangular pulses emitted by each spike $V_{\text{pre}}$. Their width is $\tau_a$ (see Sec. II B) and their intensity is controlled by a resistor $R_W$ in a conventional current-source implementation. The technical details are provided in the Appendix. Here, it is important to mention that $R_W$ controls the synaptic coupling that is denoted $W_{ij}$ between neurons i and j, in the terminology of artificial neural networks.

These $W_{ij}$ couplings are key parameters for networks that implement synaptic learning, such as through spike time-dependent plasticity [3]. In that case, one may adopt nonvolatile memristors for the $R_W$, as we shall discuss elsewhere.

The synaptic current pulse $I_{\delta}$ has a short duration, so it approximates a delta function $I_{\delta}$ that is leaky integrated by the $R_s C_s$ pair. This circuit provides the time constant $\tau_s$ and implements the differential equation in the current $I_{s1}(t)$, which is the output of the second block.

c. The third block. C implements Eq. (12) in a similar fashion. There, the current $I_{s2}(t)$ is the dynamic variable and the current $I_{s1}$ enters as the independent term. To do this we need to take $I_{s1}$ output from the second block and input it to the third one. The tricky part is to “copy the current” from one to the other, independently of the load represented by the latter. We resolve this by a standard electronic circuit, aptly called a “current mirror” (CM), and provide the technical details in the Appendix.

We note a useful feature of our CM implementation, namely, given the input, the output can be chosen either positive or negative in polarity. This feature is extremely convenient, as we can implement either excitatory or inhibitory synaptic currents. In other words, this feature allows one to choose the sign of the synaptic coupling parameters $W_{ij}$ in an artificial neural network.

This circuit block leaky integrates the input current $I_{s1}$ with timescale constant $\tau_s = R_s C_s$, similarly as the previous stage. The output of this block is the current signal $I_{s2}(t)$ and the same CM solution is adopted to input it (or its negative copy) to the last block, which is the aLIF circuit of the postsynaptic neuron.
d. The last block. D implements Eqs. (9) and (10), i.e., the aLIF neuron model [28]. This subcircuit has already been discussed in detail above. Here we just note that it corresponds to the postsynaptic neuron, so the action-potential spike emitted by this last block is denoted by $V_{post}(t)$.

## 2. Synaptic current forms

The reader may have realized that our implementation provides the freedom to adopt either one or two leaky integration synaptic blocks. If only one (yellow in Fig. 3) block is used, the single spike input $I_δ$ is transformed into $I_{s1}(t) \sim \pm \exp(-t/\tau_s)$. While, if two blocks are adopted, the resulting output to a single input $I_δ$ spike is $I_{s2}(t) \sim \pm \alpha(t/\tau_s)$.

In Fig. 4, we show an example of the temporal traces of the currents at different points of the NS-unit circuit of Fig. 3. The first one is the presynaptic action-potential spike, the output of box A in Fig. 3, which shows the

<figure><img src="image_6.png" /><figcaption>FIG. 4. Traces at the different stages of the NS unit. (a) The presynaptic neuron action-potential spike is the input to the NS unit. (b) The action potential is transformed into a $δ(t)$-like pulse with an intensity $I_δ$ (modulated by $R_W$) and a duration given by the spike emission time $\sim \tau_a$. The pulse carries a total “neuro-transmitter” charge $Q_N \approx I_δ\tau_a$. (c) First leaky-integration stage [Eq. (11)], which provides and exponential current $I_e=I_{s1}$ with timescale $\tau_s = R_sC_s$, where $Q_N$ is the time integral of the current. (d) Second leaky-integration stage [Eq. (12)], which provides an $\alpha(t)$-like synaptic current $I_α = I_{s2}$. Note that, as in the mathematical model, both $Q_N$ and $\tau_s$ remain constant (to a reasonable approximation) after each successive leaky integration. They are also in reasonably good agreement with the nominal parameter values $Q_N \approx 80$ nC [calculated from Eq. (15)] and $\tau_s = 1.6$ ms.</figcaption></figure>

triangular shape produced by the fast discharge of $C_m$ on $R_a$ with $\tau_a \sim$ milliseconds, as we described before. This presynaptic potential, through the condition $V_{pre} > \Omega \approx 0.7$ V in our implementation (see the Appendix), produces the $I_δ(t)$ current pulse, which is input to box B. It has duration $\sim \tau_a$ and intensity $I_δ \approx V_+/R_W$, where $V_+ (\approx 4$ V) is set by the circuit voltage source (see the Appendix). The sharp pulse is then leaky integrated [cf. Eq. (11)] to produce the synaptic current $I_{s1}(t)$ with the exponential form as in Eq. (11), at the output of box B. The second leaky integration produces the synaptic current $I_{s2}(t)$ with the alpha form as in Eq. (12), at the output of box C. As shown by the good quality of the fits in Fig. 4, the NS faithfully implements those functional forms and keeps constant the normalization of the integrated charge $Q_N$ to a reasonably good approximation.

## 3. Synaptic current versus conductances

There is one important point about our implementation that is worth clarifying here. Theoretical models are formulated in terms of either membrane conductances or synaptic currents (see Ch. 3 and Fig. 3.2 in Ref. [3]). However, we shall show that these two options are essentially the same.

The key point is to realize that neurons spend most of their time at and around their resting potential, $V_m \approx V_{rest} \approx -70$ mV. The potential $V_m$ slowly drifts up or down under the influence of synaptic currents. Eventually, a small upward change of about 10 mV drives the neuron to emit a short action-potential spike and quickly returns to $V_m \approx V_{rest}$. The general expression of synaptic currents in theoretical models reads [3]

$$I_{\text{syn}}(t) = -g_{\text{syn}}(t)(V_m(t) - E), \quad (13)$$

where $g_{\text{syn}}(t)$ is the synaptic conductance and $E$ is a constant potential parameter, which is chosen either higher or lower than $V_{rest}$ for excitatory or inhibitory currents, respectively [3,13]. Hence, from $V_m(t) \approx V_{rest}$ all the time except in brief firing events lasting a few milliseconds, to a good approximation, the temporal dependence of $I_{syn}$ is the same as $g_{\text{syn}}(t)$, so one can model either. This feature is also the reason why one often reads in the neuroscience literature expressions such as “to inject a conductance” [29].

## D. Synaptic current timescales: spike by spike to rate coding

An important question in neuroscience is what is the nature of the neural code. In neural networks, do neurons code through the timing between spikes or by the neuron spiking rate? Our neurosynaptic unit is well equipped to implement both of these representations. The key parameters to compare are the timescale for successive spike
<figure><img src="image_7.png" /><figcaption>FIG. 5. (a) Spike-by-spike coding. Each emitted presynaptic spike (top) produces individual α-type spikes of synaptic current (bottom). (b),(c) Rate coding. The presynaptic emitted spikes (top) have a high rate and produce a synaptic current buildup by accumulation or leaky integration (bottom). The current intensity encodes the spiking frequency rate, $I_s = I_s(f)$.</figcaption></figure>

emission, i.e., the membrane timescale $\tau_m \sim R_m C_m$, versus the synaptic integration timescale $\tau_s \sim R_s C_s$.

If $\tau_s \le \tau_m$ then the synaptic leaky integration can follow individual spike emission and spike-by-spike time coding is possible. On the other hand, if $\tau_s \ge \tau_m$, the effect of the successive spikes gets accumulated and a synaptic current $I_s$ proportional to the firing rate $f$ results. In Fig. 5 we illustrate the different coding regimes. We show the input spikes $V_{pre}$ of a presynaptic neuron acting on the synaptic unit (cf. Fig. 3) and the resulting output synaptic currents $I_{s2}$ that it generates.

To determine whether the system is in the spike-by-spike or rate-coding regime, we need to consider the activation function of the neuron model $f(I_{in})$ that is shown

<figure><img src="image_8.png" /><figcaption>FIG. 6. Response functions of neuron and synapse subcircuit units. (a) The spiking rate as a function of the input current, $f(I_{in})$, produced by the aLIF neuron subcircuit (box D in Fig. 3). Here $I_{min}$ denotes the threshold of excitation and $I_{max}$ the value where the neuron stops spiking. (b) The synaptic current $I_s(f)$ of the leaky-integrator subcircuit (boxes B and C in Fig. 3) in the spike rate regime. The slope is given by the neurotransmitter charge $Q_N$, as in Eq. (14). The dashed line shows the theoretical $I_s(f)$ for the nominal value of $Q_N$ computed with Eq. (15).</figcaption></figure>

in Fig. 6. The frequency or equivalently the interspike interval $\text{ISI} = 1/f$ depends on the external input current excitation, so it is more accurate to compare this timescale (instead of $\tau_m$) with the leaky-integration timescale $\tau_s$. Hence, per our previous discussion, we have that $1/f > \tau_s$ would correspond to spike-by-spike coding, while $1/f < \tau_s$ would correspond to the rate-coding mode. Thus, higher spiking frequencies naturally correspond to the rate-coding regime, and this is well illustrated in Fig. 5.

We may make here a side comment on the activation function $f(I_{in})$. As shown in the Fig. 6 this function also exhibits a maximal current threshold $I_{max}$, where the firing frequency suddenly decreases down to zero [28]. This is actually a biomimetic property of our neuron model. The eventual decrease of the spiking frequency at high current excitation is a phenomenon that is observed in biological neurons, where is known as “depolarization block” [30]. It is often considered relevant for the understanding of epileptic seizures [31].

It is an important question to determine what controls the intensity of the synaptic current in the rate-coding mode. We have discussed in Sec. IIC2 how the intensity of the synaptic current that is produced by a single action-potential spike is controlled by the resistor $R_W$. More specifically, we showed that for spike-by-spike coding, $R_W$ controls the total charge due to a single spike, i.e., $Q_N = \int I_s dt$. This charge per spike is conserved across the leaky-integration stages, i.e., there is no current amplification.

For the rate-coding case, the successive spikes are accumulated, hence we expect that the synaptic current will be given by $Q_N f$. A subtle point to note, however, and as seen in Fig. 5, is that the accumulation takes a certain time to reach a steady state, which is the leaky-integration time constant $\tau_s$. Therefore, after the transient integration time, we have

$$I_{se} \approx Q_N f \qquad (14)$$

for the excitatory case. This simple relation is in good agreement with our circuit measurements shown in Fig. 6, including the quantitative value of the slope $Q_N \approx 80$ nC (see Fig. 4).

We may provide the explicit dependence of $Q_N$ with the circuit model parameters. Recalling the expression for the charge contribution of each spike, we have

$$Q_N \approx I_\delta \tau_a \approx (V_+/R_W)R_a C_m. \qquad (15)$$

Finally, we should note that all the previous considerations also remain valid for the case of inhibitory synaptic currents (since we just change the polarity of the last current mirror). Hence, in that case the $I_s(f)$ response function is simply the reflection of Eq. (14), namely,

$$I_{si} \approx -Q_N f. \qquad (16)$$
### III. RESULTS

#### A. Basic spiking-neural-network motifs and neurocomputational primitives

The methodology that we described above, in terms of circuit blocks, is modular, and therefore is well adapted to implement arbitrary biomimetic SNNs. In the next sections, we shall illustrate how one can establish basic network motifs, such as a single neuron with synaptic self-excitation or self-inhibition, and then two mutually inhibiting neurons. The spiking states of these motifs realize neurocomputational functions or primitives.

The basic motifs that we shall consider are shown in Fig. 7. The simplest motif is, of course, an isolated spiking neuron. We have already discussed its behavior in Sec. II B. The first nontrivial motifs display self-feedback, which is the simplest case of a recursive SNN. The self-excitatory case provides an instantiation of an important neurocomputing primitive, called dynamical memory, short-term memory, or simply activity “bump.” The second case is a self-inhibitory neuron, which implements two other important neural primitives: spiking adaptation and bursting [3,13].

Then, we shall leverage those basic primitives to build more complex two-neuron SNNs. We shall consider two different motifs, and implement three different types of central pattern generators (CPGs) [13,32]. These are also fundamental neuronal functions that are widely studied, as they broadly appear in neuroscience [33,34] and robotics

<figure><img src="image_9.png" /><figcaption>FIG. 7. SNN motifs. (a) Single spiking neuron. (b) Single spiking neuron with excitatory feedback. This motif represents both an autapse (self-axon-dendrite connection) or an intrinsic excitation current, like the inward $Ca^{2+}$, $I_T$. (c) Single spiking neuron with inhibitory feedback. This motif represents both an autapse or an intrinsic adaptation current, like the hyperpolarizing $I_{AHP}$ that is also mediated by $Ca^{2+}$ (see Sec. 10.3 of Ref. [13]). (d) CPG primary-subordinate motif with a burster neuron (primary) that inhibits a regular spiking neuron (subordinate). (e) CPG of mutually inhibiting neurons with self-adaptation. The same motif has two cases: two bursters or two frequency-adaptation neurons, which give qualitative different traces.</figcaption></figure>

[35,36]. While in the present introductory work, we shall stop at two-neuron motifs, it should be clear that the present methodology is general and permits building larger SNNs with arbitrary architecture.

#### B. Self-excitation: Activity bump and dynamical memory

In this section, we begin with the simplest motif: a self-excitatory neuron (see Fig. 7). This simple motif brings a significant surprise: contrary to the common wisdom, it is sufficient to implement a basic neurocomputational primitive, namely, a dynamical memory [3]. Indeed, it is expected for a single neuron with a self-excitatory feedback to show a runaway rate, and that a population of neurons is needed to achieve a stable self-excitation [3]. Here, we shall demonstrate that this is not the case by explicitly exhibiting the dynamical memory in hardware and also its stability against strong perturbations.

The reverberating self-sustained state is of central importance in neuroscience, where it is called short-term memory (cf. Ch. 6 in Ref. [13]) or working memory [37]. The excitation “bump” neurocomputational function occurs in models of great current interest, such as place and grid cells [38] and in head direction systems [39]. It also appears in models of the oculomotor system, which allows the gaze to be fixed in a given direction [40].

In a biological neuron, the self-excitatory current may represent the intrinsic transient inward $Ca^{2+}$ ionic current, known as $I_T$ (see Ch. 10 in Ref. [13]). If, alternatively, the self-synaptic connection is extrinsic, such as due to the neuron's own axon connecting the dendrites, it is called an autapse, which is a topic of current interest [41-44]. This may be relevant to epilepsy, where it was reported that, in epileptic humans, about 30% of cortical neurons of layer V form autapses [45]. Moreover, very recently, it was observed that the ictal source point of epileptic seizures occurs in layers IV, V, and VI [46].

For a self-excitatory synaptic current, the system of equations (9)–(12) of our model become

$$C_m \frac{dV_m}{dt} = -\frac{V}{R_m[S] + R_a} + I_{sse} + I_0, \quad (17)$$

$$V_{ap} = \frac{V_m}{R_m[S] + R_a} R_a, \quad (18)$$

$$\tau_s \frac{dI_{s1}}{dt} = -I_{s1} + I_{\delta} H(V_{ap} - \Omega), \quad (19)$$

$$\tau_s \frac{dI_{sse}}{dt} = -I_{sse} + I_{s1}. \quad (20)$$

Here, for the self-excitation, the spike potentials become $V_{pre} = V_{post} = V_{ap}$. Additionally, to indicate the synaptic self-excitatory current, $I_{s2}$ is renamed $I_{sse}$. The system of equations above corresponds to the neurocomputing circuit shown in Fig. 8. The key feature of the activity bump
<figure><img src="image_10.png" /><figcaption>FIG. 8. Self-synaptic current circuit. Here $I_0$ is the external current input; and $V_{ap}$ is the output train of action-potential spikes. The output spikes are feedback to the neuron through the self-synaptic current $I_{sse}$.</figcaption></figure>

is that the neuron can reach a spiking state, and remain self-excited, despite the input external current $I_0$ being sub-threshold, namely $I_0 < I_{min}$. As we shall discuss below, the self-excitation state is maintained by the buildup of the self-synaptic current $I_{sse}$.

To start the bump of activity, it is necessary to apply a short excitatory external current pulse $I_\Delta$ on top of the constant subthreshold current $I_0$, as shown in Fig. 9. Since $I_0$ is subthreshold, we observe that there is no activity during the first 2 s. Then an excitatory pulse is applied, that fulfills two requirements: (i) it has to drive the neuron to the excitatory regime above threshold, i.e., $I_{min} < I_0 + I_\Delta < I_{max}$; and (ii) its duration $\tau_\Delta$ has to allow for the buildup of the self-synaptic current $I_{sse}$, i.e., $\tau_\Delta > \tau_s$.

We observe in Fig. 9 that the spiking starts as soon as $I_\Delta$ begins, and within the next second ~ $\tau_s$ the $I_{sse}$ starts to build up and the frequency rate increases. Within the next second, by the time $I_{sse}$ is terminated, the self-feedback is well in place. The persistent state is realized as $I_0 + I_{sse} > I_{min}$. The spike rate of the bump of activity stabilizes at the self-consistent rate $f^*$ as we discuss below.

We can recall the two response functions that characterize the neurosynaptic unit (see Fig. 6) and use them to formulate the bump state as a self-consistent problem. On the one hand, we have $f(I_{in})$, the neuron activation function; on the other, the synaptic current $I_s(f)$ that results from the synaptic circuit, given an input spike rate $f$ from a neuron. We should also note from Fig. 9 that $\tau_s > 1/f$, so we are indeed in the rate-coding regime.

We would now like to analyze the dynamical memory state under a different light. We shall show that it can be cast as a self-consistent problem, which can be solved by a geometrical construction. The same method will be used in the next section to consider the case of self-inhibition, i.e., adaptation.

In the self-excitatory motif, the synaptic current is feedback to the neuron, hence $I_s = I_{sse}$ (see Fig. 6). Therefore,

<figure><img src="image_11.png" /><figcaption>FIG. 9. (a) Top: Spike trace of the excitation bump. The blue line schematically denotes the applied external subthreshold current $I_0$ with the short excitation $I_\Delta$ that starts the bump, and inhibition that terminates it. Bottom: The instantaneous spike-rate emission $f$ of the neuron (black line); and the self-synaptic current $I_{sse}$ buildup during excitation and relaxation during inhibition (red line). The state has an external continuous subthreshold $I_0 = 44.7$ µA. (b) Geometrical construction to determine the self-consistent frequency $f^*$. Insets: $f^*$ as functions of $R_w$ and $I_0$. The parameters are $\tau_s = 0.47$ s, $I_0 = 44.7$ µA, $+I_\Delta = 22.9$ µA, $-I_\Delta = -24.4$ µA. For clarity, the number of spikes is reduced to one-tenth of the originally measured number.</figcaption></figure>

we can cast the bump as a self-consistent problem: (i) the total input current is $I_{in} = I_0 + I_{sse}[f]$, where $I_0$ denotes a constant subthreshold external current and $I_{sse}$ is the self-excitatory contribution that depends on the spiking frequency; and (ii) the total $I_{in}$ produces a spike rate $f(I_{in})$, which itself produces the feedback current $I_{sse}$. Therefore, the problem of simultaneously fulfilling (i) and (ii) is to find the self-consistent frequency $f^*$ such that $f^* = f[I_0 + I_{sse}(f^*)]$.

The solution to this problem can be obtained by a geometrical construction. It corresponds to the crossing point $(I_{in}^*, f^*)$ between the two response functions measured on the NS unit (see Fig. 6). From the first response, we call
$f_1 = f_1(I_{in})$, which is shown by the blue line in Fig. 9(b). From the second response, we invert the variables to get $f_2 = f_2(I_s) = f_2(I_{sse})$, shown in dashed yellow line in Fig. 9(b) (and compare to Fig. 6). From (i) above, we have $I_{in} = I_0 + I_{sse}$; then we get $f_2 = f_2(I_{in} - I_0)$, which is the yellow solid line in Fig. 9(b). The geometrical solution is obtained by requiring $f_1 = f_2$, which gives the crossing point $(I_{in}^*, f^*)$, as shown as a red cross in Fig. 9(b).

We can further explore how the self-consistent firing rate $f^*$ systematically depends on two main parameters of the model, as we show in Fig. 9. One is the synaptic weight $R_W$, which controls the intensity of the feedback current. This intensity is controlled by the “neurotransmitter” charge $Q_N$, which is modulated by the synaptic resistor $R_W$ (see Fig. 4 and Sec. II D). The other parameter is the external current intensity $I_0$, which is always kept subthreshold for the activity bump. One could in principle set $I_0 = 0$ and still realize a bump. The requirement is a strong enough feedback, such that $I_{sse} > I_{min}$. The continuous dependence of the $f^*$ on the model parameters is shown in the insets of Fig. 9(b).

One of the long-standing issues in theoretical models of this type of dynamical attractor concerns their stability with respect to perturbations, since it is based on a positive feedback loop [19,40]. Therefore, it is important to note a conceptual point in this regard. The methodology that we adopt in the present work is based on the *physical* implementation of a theoretical model. Thus, the stability of the dynamical attractor model is directly demonstrated *by construction*. Namely, if it works, therefore, it must be stable. Nevertheless, we may further demonstrate this feature by explicitly applying strong external perturbations to the activity bump state. This is shown in Fig. 10. Despite

<figure><img src="image_12.png" /><figcaption>FIG. 10. Stability of the bump attractor state. Top: An excitatory current pulse of intensity $+I_Δ = 22.9$ µA (blue line) is applied during 1.2 s, producing an increase in the firing rate, which rapidly relaxes after the perturbation stops. Bottom: An inhibitory perturbation $-I_Δ = -24.4$ µA is applied during 0.8 s, producing a halt on the spiking. The bump recovers rapidly after the perturbation stops. Inset: Evolution of the synaptic currents $I_{s1}$ and $I_{s2} = I_{sse}$ during the perturbation. The parameters are $τ_s = 0.47$ s and $I_0 = 44.7$ µA. The number of spikes is reduced to one-fifth for the sake of clear display.</figcaption></figure>

these strong perturbations, even halting the spikes in the inhibitory case, the dynamical attractor returns promptly to its stable state, as shown in the inset. A key feature for the stability, however, is that the perturbations should not last much longer than the synaptic timescale $τ_s$. Hence, we may say that this state can keep a short-term memory encoded by $f^*$ during that characteristic timescale.

Finally, from these stability considerations, it follows that, to stop the active bump state, one needs to apply an inhibitory current pulse $-I_Δ$ that is sufficiently long and strong. More precisely: (i) its magnitude has to be such that the total input current $I_{in} = I_0 + I_{sse} - I_Δ < I_{min}$, i.e., falls below the onset of excitability, so stops the firing and allows the feedback to relax; and (ii) its duration, therefore, should be longer than $τ_s$, such that the self-synaptic current $I_{sse}$ relaxes sufficiently, as shown in Fig. 10, such that the condition (i) is fulfilled.

### C. Self-inhibition: Adaptation and bursting

In this section, we turn to another minimal motif, namely that of neuron self-inhibition (see Fig. 7). As we shall see, this motif produces two important spiking modes: adaptation and bursting. These modes will be the basis for the multineuron motif for central pattern generators, which we shall describe later in the next section.

#### 1. Adaptation

Adaptation is a relevant function in neuroscience. An example of this type of self-inhibition is the after-hyperpolarization current $I_{AHP}$, which is mediated by $Ca^{2+}$ (see Ch. 10 in Ref. [13]). Similarly, as mentioned in the previous section, this motif may also be realized by an external connection between the neuron's own axon and its dendrites. This case is called an inhibitory autapse, which is now known to be a common feature of some neocortical and parvalbumin (PV) neurons, with the function of regulating the rate of spike emission [42,45,47,48]. Self-inhibition is also a main ingredient of two of the most popular schematic theoretical models, Izhikevich and AdEx [3,17]. In both models, the self-inhibition is described by a second dynamical equation for the adaptation or recovery variable.

Similarly as done in the previous section, we shall also show how the spiking states of this motif can be cast as a self-consistent problem. The geometrical solution of the self-consistent problem brings valuable insights, specially for the understanding of the stability of the firing states and the nature of the bursting one. Rather surprisingly, despite a very large literature on the Izhikevich and AdEx models, we are not aware of any similar treatment to the one we introduce here [3].

We first begin by exploring the dynamical behavior of this important motif by means of our hardware implementation. In this case, the system of equations (9)–(12) of our
theoretical model becomes:

$$C_m \frac{dV}{dt} = - \frac{V}{R_m[S] + R_a} - I_{\text{ssi}} + I_0, \quad (21)$$

$$V_{\text{pre}} = \frac{V}{R_m[S] + R_a} R_a, \quad (22)$$

$$\tau_s \frac{dI_{s1}}{dt} = -I_{s1} + I_{\delta}H(V_{\text{ap}} - \Omega), \quad (23)$$

$$\tau_s \frac{dI_{\text{ssi}}}{dt} = -I_{\text{ssi}} + I_{s1}. \quad (24)$$

Here we now need to excite the aLIF neuron with a suprathreshold external current $I_0$, unlike in the previous section. Also, we rename the synaptic current $I_{s2}$ as $-I_{\text{ssi}}$ to denote the synaptic self-inhibition. The circuit configuration is identical to that shown in Fig. 8 for the self-excitation, with the sole difference that the mirrored current $I_{s2}$ is now of inverted polarity.

While it is possible to have some degree of adaptation in the spike-by-spike regime, when $\tau_m > \tau_s$, the most natural, interesting, and biologically relevant case is the rate-coding regime, i.e., for $\tau_m \sim 1/f \ll \tau_s$. Hence, we focus here on the latter. In Fig. 11 we show the behavior of the aLIF with a self-adapting current. The adaptation results from the buildup of the inhibitory current intensity $I_{\text{ssi}}$ on the synaptic subcircuit, which occurs within the synaptic timescale $\tau_s$. The total current exciting the neuron gets reduced to $I_{\text{in}} = I_0 - I_{\text{ssi}}$, which leads to the reduction of the firing rate with respect to the original one, $f(I_0)$. This phenomenon of spiking rate reduction is called frequency adaptation.

As we did before, we may also consider this state as a self-consistent problem, where we need to find the spike rate $f^*$. Now, the self-consistent condition is given by $f^* = f^*(I_{\text{in}}) = f^*[I_0 - I_{\text{ssi}}(f^*)]$. As we discussed before, in Sec. II D, Eq. (16), the response function $I_{\text{ssi}}(f) \approx -Q_N f$, i.e., identical to the response function $I_{\text{sse}}(f)$ with a change of sign. So the inhibitory current represents a discharge of $C_m$.

In Fig. 11 we show the geometrical construction for the case of a negative self-synaptic current, which follows the same steps outlined before. Similarly, as in the case of self-excitation, the self-consistent frequency can be controlled by modulating the intensity of the feedback current. This can be achieved by varying the model parameters, such as $R_W$ and $I_0$, as shown in the insets of Fig. 11.

Also similarly as done before for the bump, the stability of the adapted self-consistent state can also be tested by applying perturbations. Our results are shown in Fig. 12, where we demonstrate very good stability.

Before leaving this section we should make a final point. As already mentioned before, it is one's choice to adopt one or two leaky integration blocks. This choice can be made in either the excitatory or the present inhibitory case.

<figure><img src="image_13.png" /><figcaption>FIG. 11. Adaptation traces. (a) Top: Spiking trace showing the progressive reduction of the frequency. The blue line indicates the external current input $I_0 = 92.8$ µA. Bottom: The instantaneous spike rate $f$. The red line shows that the self-synaptic adaptation current $I_{\text{ssi}}$ builds up over a timescale $\tau_s \approx 1$ s; and $f^*$ denotes the self-consistent asymptotic frequency. (b) Geometrical construction to determine the self-consistent frequency $f^*$. Insets: $f^*$ as functions of $R_W$ and $I_0$. The number of spikes is reduced to one-twentieth for the sake of clear display.</figcaption></figure>

Below, we explicitly write the model equations, where we adopt the simplified notation of Eq. (1), for the sake of clarity, and adapt them for the case of one inhibitory leaky integration stage. The system reads

$$C_m \frac{dV_m}{dt} = f(V_m) - I_{\text{ssi}} + I_0, \quad (25)$$

$$\tau_s \frac{dI_{\text{ssi}}}{dt} = -I_{\text{ssi}} + I_{\delta}H(V_{\text{ap}} - \Omega). \quad (26)$$

The reader may verify that these equations now share the same form as the Izhikevich and AdEx models [3,17]. The $I_\delta$ term represents the spikes due to the action potentials $V_{\text{ap}}(t)$, which are analogous to the potential spikes of fast variable $V(t)$ in the theoretical models. In all cases, these
<figure><img src="image_14.png" /><figcaption>FIG. 12. Stability of the frequency-adaptation state. Top: An excitatory perturbation current pulse of intensity +I$_{Δ}$ = 22.6 µA (blue line) is applied during a short period. This produces an increase in the firing rate, which rapidly relaxes after the perturbation stops. Bottom: An inhibitory perturbation −I$_{Δ}$ = −24 µA (red line) is applied, producing a halt on the spiking. The adapted state recovers rapidly after the perturbation stops. Inset: Evolution of the synaptic currents I$_{s1}$ and I$_{s2}$ = I$_{ssi}$ during the perturbation. The parameters are τ$_{s}$ = 1 s and I$_{0}$ = 92.8 µA. The number of spikes is reduced to one-twentieth for clarity.</figcaption></figure>

spikes are identically leaky integrated by the recovery variable.

## 2. Bursting

We have adopted the intuitive perspective of considering states as the solution of self-consistent problems. Therefore, it is interesting to pose the following question: What would happen if we increase the intensity of the adaptation, so that the resulting total input current to the neuron falls below the excitability threshold? Namely, I$_{ssi}$ is such that I$_{in}$ = I$_{0}$ − I$_{ssi}$ (f) < I$_{min}$, so that the strong self-inhibition would stop the firing, hence the self-consistent solution is no longer possible.

Let us analyze what would be the expected behavior. The self-inhibition takes a time τ$_{s}$ to build up and stop the firing. However, this quiescence provokes the relaxation decay of the feedback current, within a timescale τ$_{s}$ . Recalling that I$_{0}$ > I$_{min}$ , the total input I$_{in}$ = I$_{0}$ − I$_{ssi}$ will grow as I$_{ssi}$ decreases and, eventually, the neuron will get re-excited and start spiking again. Thus, this dynamical bistability points to periods of spike emission alternating with quiescent ones. Such a state is called bursting, which is another basic neurocomputing primitive of the self-adaptation motif.

The bursting state is implemented by the NS unit as we show in Fig. 13. From the qualitative discussion, it should be clear that the bursting state is controlled by the two timescales τ$_{m}$ and τ$_{s}$ . The latter controls the alternation period or burst cycle, and the former the interspike intervals during each burst. Hence, we have τ$_{m}$ ≪ τ$_{s}$ so the neuron is in the rate-coding regime. Moreover, from the previous discussion, one should expect that the intraburst spiking frequency has to evolve around the excitability

<figure><img src="image_15.png" /><figcaption>FIG. 13. Burster or pacemaker state. Top: Spike traces showing bursts of activity separated by quiescent states. Middle: The total input current I in is oscillating narrowly around the threshold I min . Bottom: The self-inhibitory current I ssi resembling a harmonic oscillation with period ≈ τ s = 1 s, as discussed in the text.</figcaption></figure>

threshold I$_{min}$ . As shown in Fig. 13 this is indeed the case, with I$_{ssi}$ oscillating around the threshold I$_{min}$ within a fraction of a µA.

We may gain further and deeper analytic understanding of those small oscillations. This is revealed by a simple analysis of the model equations, which to our knowledge has not been done before. We consider the model equations and take the time derivative of Eq. (24) and replace it into Eq. (23) to obtain

$$
\tau_s^2 \frac{d^2 I_{\text{ssi}}}{dt^2} + 2\tau_s \frac{dI_{\text{ssi}}}{dt} + I_{\text{ssi}} = I_{\delta}[f(t)]. \quad (27)
$$

This differential equation is the analogue to the familiar driven damped harmonic oscillator model [49] with resonant frequency 1 /τ$_{s}$ . The “friction” or damping is the second term and the “external driving force” is I$_{δ}$ [f (t)]. From the behavior of this well-known system, we may expect two different states: one where the damping dominates, and another where the oscillations persist. Indeed, as shown in Fig. 13, the synaptic current I$_{ssi}$ (t) exhibits a beautiful oscillatory behavior with a period close to τ$_{s}$ = 1 s, as predicted by Eq. (27).

However, our driven damped oscillator has a twist, since the driving force is not arbitrary but depends on the neuron’s spiking trace through the feedback loop. So we can discuss the nature of the oscillation in more detail. Let us assume that the system is bursting and show that this assumption is consistent. The current I$_{δ}$ [f (t)] consists of trains of fast spikes at frequency f separated by silent periods, thus having a rhythmic time dependence with period τ$_{s}$ . The trains of current spikes are twice leaky integrated with timescale τ$_{s}$ , therefore, I$_{ssi}$ (t) builds up within the
same timescale as the period. Similarly, during the quiescent phase, the inhibitory current relaxes on the same timescale. Then, the periodic buildup and relaxation of $I_{ssi}(t)$ modulates the firing behavior of $f(t)$, hence of $I_δ$, and the bursting state is consistent. The point to make is that the “driving force” $I_δ(t)$ is in resonance, as the trains of spikes are emitted at a frequency $1/τ_s$, which is the natural frequency of the oscillator equation (27). Within the same line of reasoning, the simple adaptation can be understood as an overdamped oscillation.

## IV. CENTRAL PATTERN GENERATORS

We are now ready to start building general functional spiking neural networks. We shall begin by considering the simplest multineuron motif, namely two neurons coupled by dynamical synaptic currents [20]. We shall focus on implementing an SNN whose neurons emit periodically alternating bursts. Such a two-neuron network is called a central pattern generator and has a very important role in neuroscience [32]. This type of spiking activity is fundamental for motor systems with periodic dynamics, such as heartbeats, walking, eating, etc. Moreover, CPGs are also a key component in robotics, self-propelled systems, and biomedical devices. Therefore, our simple and systematic methodology may open an interesting avenue for applications. We want to emphasize here that the methodology of combining motifs is general. One can go on combining motifs in a modular manner to achieve more complex spiking behaviors, such as, for instance, the classic problems of gaits and swimming [35]. A technical side remark is to realize here that the formulation of our circuit models is done in terms of currents and implemented by current sources and current mirrors. Unlike voltages, input current sources can be added algebraically.

The CPGs involve an external excitation that induces the basic tonic spiking of each neuron, which is then sculpted by inhibitory synaptic currents acting on different timescales. As this is a vast field, we shall adopt the review of Marder and Bucher [32] as guidance to select a few paradigmatic systems to build. Nevertheless, in doing that, we shall also obtain insights provided by the present hardware-based approach.

We shall consider three different basic CPGs [32] that correspond to two different two-neuron motifs, those in Figs. 7(d) and 7(e). One has primary-subordinate or feed-forward inhibition architecture, while the other has recurrent connections of mutual inhibition. In the model circuits here, we adopt all synaptic currents, intra- and interneuron, as inhibitory $α(t)$-type functions.

### 1. Primary-subordinate CPG

We start with the simplest case, the motif in Fig. 7(d). Such CPG has a bursting neuron (N1), acting as a

<figure><img src="image_16.png" /><figcaption>FIG. 14. CPG primary-subordinate. Top: Uncoupled case, in which N1 is a pacemaking burster and N2 is a regular tonic spiking neuron. Bottom: Coupled case, in which the N2 subordinate is inhibited by the activity of N1, creating a complementary bursting pattern. Parameters are $τ_{ssi} = 1$ s for the pacemaking and $τ_{spi} = 1$ ms for the fast inhibitory projection.</figcaption></figure>

pacemaker, that sculpts the firing rate of a regular spiking neuron (N2) by projecting inhibition. The pacemaking neuron is realized by a bursting neuron, as discussed in the previous section. The subordinate neuron shows regular spiking when uncoupled, as shown in Fig. 14. When the inhibitory projection is introduced, the subordinate neuron develops a burst spiking pattern that is the inverse of its primary’s (see Fig. 14). Hence, the two-neuron system displays a sequence of alternating bursts which is characteristic of CPGs.

A requirement to achieve this state is that the interneuron synaptic projection from the primary has to be strong enough to silence the subordinate, i.e., to drive the subordinate under its excitation threshold. So for the currents acting on N2, we have

$$I_{in} = I_0 - I_{spi} < I_{min} \quad (28)$$

where $I_{spi}$ denotes the synaptic projected inhibitory of N1 on N2. The intensity of the projection can be easily controlled by the resistance $R_W$ of the $I_{spi}$ synaptic circuit.

### 2. Mutual-inhibitory CPG

We turn now to the second motif, with mutual inhibition. In this case, we shall consider two qualitatively different possibilities: two coupled bursters and two coupled adaptive neurons. For simplicity, we shall consider symmetric systems. However, it is important to acknowledge that inherent variability exists in all hardware as a result of manufacturing tolerances. Nevertheless, as with biological neural neurons that are also nonidentical, the emergent function has to be robust to that type of variability.

Two mutually coupled pacemakers present a more straightforward case; therefore, we begin our discussion with this scenario. We consider two neurons that, if decoupled, both produce the same bursting (cf. Sec. III C 2), up to the intrinsic variability just mentioned. Therefore, the actual observed behavior is of two similar, but not identical, periodic burstings, as shown in Fig. 15. By introducing
<figure><img src="image_17.png" /><figcaption>FIG. 15. CPG mutual inhibitory pacemaker. (a) Uncoupled case: N1 and N2 are pacemaking bursters with nominally equal, but in practice similar, frequencies. (b) Coupled case: N1 and N2 settle into a common antisynchronous bursting pattern. (c) Synaptic currents of N1: the self-inhibitory $I_{ssi}$ with $\tau_{ssi} = 0.47$ s (red) and the N1 to N2 projected inhibition $I_{spi}$ with $\tau_{spi} = 0.1$ s (blue). The green line marks the beginning of an N1 burst. (d) The same as panel (c) but for the synaptic currents of N2. For clarity, the number of depicted spikes is reduced by one-fifth of the originally measured number.</figcaption></figure>

coupling in the form of projection of mutual inhibition [see the motif in Fig. 7(e)], we create anticorrelation in the respective firing patterns. As we can see in Fig. 15(b), the state of the two-neuron network finds a compromise where they burst in antiphase. Therefore, they both lock into a single common frequency. This emerging frequency is close to the uncoupled ones, but relatively reduced due to the global effect of additional inhibition.

In this case, it is interesting to observe the synaptic currents, which we show in Fig. 15. They bring interesting insights. Following the green dashed line we can correlate the burst of each neuron with the respective effects on their synaptic currents. The first thing to realize is that the burst induces $\alpha(t)$-type currents, which can be clearly appreciated. The maximum of $\alpha(t)$ indicates the respective time constants, which we denote $\tau_{ssi} = 0.47$ s and $\tau_{spi} = 0.1$ s for synaptic self-inhibitory (red) and synaptic projected-inhibitory (blue).

If we focus on the synaptic currents of N1, shown in Fig. 15(c), we can see that the burst starts when the self-inhibitory $I_{ssi}$ decreases to almost its lowest intensity. Conversely, the burst is terminated shortly after due to the buildup of $I_{ssi}$, i.e., occurring within a timescale somewhat shorter than $\tau_{ssi}$. On the other hand, we observe that the antisynchrony is due to the strong inhibitory projection, which has a relatively faster timescale $\tau_{spi}$ and prevents the superposition of bursts. Similar considerations can be made for the other neuron's currents. We can also observe that the different initial traces in the two decoupled neurons [Fig. 15(a)] also reflect some differences in the respective

produced currents. Nevertheless, the two neurons robustly lock to a common and stable pattern of alternated bursts.

It would be interesting, but beyond the scope of the present study, to perform a full exploration of the different emergent states of the two coupled bursters. Nevertheless, we can now qualitatively discuss some relevant aspects regarding the relative timescales.

In the case discussed above, we note that the duration both of the bursts and of the projected inhibition $\tau_{spi}$ are relatively short, leading to the antisynchrony already discussed. However, one may choose to increase $\tau_{spi}$, so that it is significantly longer than the bursts. Then the projection from, say, N1 to N2 will be delayed compared with the burst in N1. This projected inhibition would make the burst of N2 in antisynchrony less likely. The result is that the emergent state would be that of synchrony between N1 and N2. This situation is qualitatively similar to that analyzed by van Vreeswijk et al. [20] for the case of two simple spiking neurons coupled by a retarded mutual inhibition.

We may finally turn to another, more subtle, form of CPG, which shares the same motif as in Fig. 7(e), but is found in a different regime. This CPG is realized with two neurons that exhibit spike-frequency adaptation, i.e., are tonic spiking when uncoupled (cf. Sec. III C 1). This type of CPG model is a classic topic of neuroscience, as it finds a biological realization in the heartbeat neural activity of the leech [33,34,50].

The behavior of this mutually adaptive CPG is shown in Fig. 16. It also shows alternate burst emission, but with a different mechanism. We note that, for illustrative purposes, we have maintained the synaptic time constants the same as in the previous case. The present mechanism is achieved by mainly driving N1 and N2 to the adaptation regime, reducing the intensity of the self-inhibition (see Sec. III C 1).

As seen in Fig. 16 the CPG traces show qualitative differences with respect to the previous case. One is that the burst shows some overlap. Also in contrast to the previous case, we can now observe a significant frequency modulation, where $f(t)$ decreases at the onset and the termination of the burst. This modulation is correlated to the total input current $I_{in}(t)$, which shows a similar behavior. Note that it crosses the excitability threshold $I_{min}$ as it produces the bursts.

To appreciate the mechanism, it is essential to look at the synaptic currents generated by each neuron. They are shown in Figs. 16(e) and 16(f). As before, in red is the self-inhibition $I_{ssi}$, and in blue is the projected inhibition to the other neuron $I_{spi}$. The green dotted lines are a guide to indicate the beginning and end of the burst of N1. The analysis for N2 would be the same.

We focus on the start of the burst of N1. We observe that its self-inhibition $I_{ssi}$ is decreasing towards its minimum [red line in Fig. 16(e)], similarly as in the CPG case discussed before. However, unlike that case, now N2 is firing
<figure><img src="image_18.png" /><figcaption>FIG. 16. CPG mutual inhibitory adaption. (a) Uncoupled case: N1 and N2 are adaptive neurons with tonic spiking. (b) Coupled case: N1 and N2 settle into a common antisynchronous bursting pattern. (c) The instantaneous frequency rate of N1. The green lines denote the beginning and end of a burst. The horizontal dashed-dotted line indicates the adapted frequency $f^*$ when the neurons are decoupled. (d) Total input current $I_{in}(t)$ to N1. The beginning and end of the burst coincide with the crossing of the excitation threshold $I_{min}$. (e) Synaptic currents of N1: the self-inhibitory $I_{ssi}$ with $\tau_{ssi} = 0.47$ s (red) and the N1 to N2 projected inhibition $I_{spi}$ with $\tau_{spi} = 0.1$ s (blue). (f) The same as panel (e) but for the synaptic currents of N2. For clarity, the number of depicted spikes is reduced by one-fifth of its original measured number.</figcaption></figure>

and intensively projecting its inhibition $I_{spi}$ on N1 [blue line in Fig. 16(f)]. In consequence, the reason why N1 starts to fire is only due to the decrease of its self-inhibition. We may say that N1 recovers and escapes from the grip of N2. This is one of the four typical CPG mechanisms, and is indeed called “intrinsic escape” [23,51].

We may also look at the end of the burst phase. There we observe that both $I_{ssi}$ and the received projection $I_{spi}$ are on the rise, so they both contribute to inhibiting N1. Moreover, the projection of N2 is *necessary* to increase the inhibition of N1 and render it quiescent. This is because N1 is now an adaptation neuron and not a pacemaker as in the previous case.

This last observation is relevant because it explains an evident qualitative difference between the patterns of the two modes of CPG that we considered. By comparing the spiking traces of Figs. 15 and 16, we observe that, in the first case, the bursts do not overlap in time, while in the second they do. The reason is that, in the first case, the ability to burst is an intrinsic feature of each neuron, while the inhibitory projection solely creates an anticorrelation

between them and lengthens the period. In contrast, in the second case, the neurons do not have the pacemaking ability to quiet themselves. Thus the overlap of activity is a *necessity*. This qualitative difference is a robust feature that may help in the classification and understanding of CPGs in biological systems. For instance, the heartbeat of the leech is a classic CPG animal model and is a clear example of overlapping bursts [50]. On the other hand, the CPG of the stomatogastric ganglion of crustaceans, which is another paradigmatic animal model, is a clear example of pacemakers with anticorrelations [32].

## V. DISCUSSION: EXTENSIONS TO OXIDE MEMRISTORS AND TO VLSI

### A. A template for memristive material research

In this work, we have introduced a circuit implementation of a general-purpose, theoretically based, neurosynaptic model using solely conventional electronic components. This has the evident advantage that the methodology can be easily reproduced at a very affordable cost and wide availability.

Nevertheless, our methodology also opens an exciting perspective for researchers in material science working on memristive systems. In fact, our circuit can be adopted as a template for testing memristive devices and building small to midscale functional networks. More specifically, there are two places where memristors may find their place. The first one is in the neuron stage. There, the memristor device (thyristor plus resistor) can be replaced by a *volatile* oxide memristor, such as $VO_2$, $V_2O_3$, $V_3O_5$, $NdNiO_3$, $SmNiO_3$, $NbO_2$, etc. [26,52–56].

The synaptic resistance $R_W$ is the second place where an oxide memristor, which should now be a *nonvolatile* memristive material, can replace a conventional resistor in our NS-unit circuit. The $R_W$ resistance controls the intensity of the synaptic coupling $W_{ij}$ between two neuron units i and j, in the language of artificial intelligence systems. Nonvolatile oxide memristors such as $TiO_2$, $Ta_2O_5$, $HfO_2$, etc., are popular choices for synaptic cross-bar neuronal networks [26,57]. The plasticity in the resistance modulation of oxide memristors may further endow the network with learning ability [58].

### B. Positioning and comparison to other neuromorphic hardware

A lot of overlap and confusion exist in the field of hardware-based neuromorphic systems. So we may first distinguish our approach by enumerating what our system is not: (i) It is not based on transistors in the subthreshold regime, as CMOS chips are. Thus, it is much less affected by the variability mismatch issue. (ii) It is not based on digital hardware, which is not well adapted for the continuous nonlinear dynamics of excitability characteristic of spiking
systems. (iii) It is not based on memristive oxide systems, which are still under research in the realm of material science. (iv) It is not based on software simulations on neuromorphic digital chips, such as field-programmable gate arrays (FPGAs), Loihi, etc., which present the issue of matching biological timescales [34].

On the other hand, we should also enumerate what our system is: (i) It is a novel hardware implementation of a textbook set of equations [13] of a general spiking neurosynaptic model. (ii) It is implemented with elemental electronic circuitry (see the Appendix). (iii) It is implemented in a clear, direct, and quantitative manner, where every model parameter has a hardware correlate enabling full tunability (see Sec. II C). (iv) Its electronic components are widely available and economic, so our system can be immediately duplicated and the data reproduced (see the Appendix). (v) It provides a direct hardware implementation with continuously tuned and biocompatible timescales (see Sec. IV), so it is ready to be deployed for applications in robotics, medical devices, control systems, etc.

We would now like to comment on the size of the SNNs that we may reach with the present methodology. The limitations do not come from fan-in and fan-out, which are essentially unconstrained. This is because we formulated all our models in terms of currents which can be added at the input of neurons. Thus the limitations are due mainly to the practical size of the systems that one may want to implement. We can estimate this on the order of hundreds of units, or maybe more for simple architectures. This should be sufficient, not only for the applications mentioned above, but also for research in theoretical neuroscience. In fact, it is now increasingly recognized that the dynamics of biological neural systems can often be described in low-dimensional manifolds.

While the present approach is qualitatively different than those based on VLSI, we should mention that it can be ported to VLSI, if it is so desired. In fact, the implementation of the memristive neuron unit has already been reported [59], where the function of the thyristor was emulated by two transistors. On the other hand, the synaptic unit is based on current mirrors, which are standard to VLSI. A single stage of a leaky-integration synapse has already been implemented some time ago [12,60]. So, the present improvement based on a two-stage integration that implements Rall's function should not present any significant impediment.

Nevertheless, it is important to emphasize that two main features would still distinguish the present methodology from previous VLSI approaches. First, and more important, our neuron circuit design would likely suffer less from device mismatch, since it is based not on subthreshold behavior but on the memristive effect [59]. Second, there is the higher simplicity of our circuit designs.

Another important aspect that needs to be commented upon is the power efficiency. While the present work is not

focused on engineering aspects, this feature is relevant for eventual practical applications. The main qualitative difference from current VLSI implementations is, of course, the size of the components. One may naively assume that the miniaturization of VLSI may have a dramatic impact on power consumption. However, this is not the case. The energy cost per spike can be estimated as proportional to $CV_{th}^2$, where C is the capacitance of the membrane capacitor and $V_{th}$ is the threshold voltage for spike emission. In VLSI, the miniaturization implies small values of C, but it also implies high spike frequency rates. In fact, the energy per spike is proportional to C, but the frequency is inversely proportional to it. Therefore, the electrical power (energy of a spike times the spiking frequency) becomes independent of C, hence miniaturization has no effect. On the other hand, our bottom-up method can use a smaller number of spiking units compared to VLSI, which needs to average populations of neurons to mitigate the mismatch [15]. Thus, achieving the same function with a dramatically smaller number of units does imply lower power requirements for our approach.

To compare the power efficiency with digital computers is less easy to establish. A practical possibility to address this question would be to compare the implementation of our CPGs with that achieved using FPGAs [34]. This is an interesting issue to address, which lies beyond the scope of the present work.

## VI. CONCLUSIONS

In the present work, we have introduced a novel neurosynaptic device that is the building block of a methodological platform to implement general spiking neural networks with a bottom-up approach. More specifically, we introduced a synaptic unit with biomimetic dynamics that, combined with our previously introduced spiking neuron [27,28,61], completes a platform to build arbitrary SNNs in hardware.

We illustrated the advantages of the approach by considering several basic motifs of theoretical neuroscience, whose dynamics has already produced results and insights.

Firstly, we found that a dynamical memory can be implemented with just one recursive spiking neuron, in contrast to the common wisdom that a population is required. The stability of the memory state is validated explicitly by the methodology that is based on measurements done on hardware implementation.

Secondly, we introduced a model of self-adaptation, similar to Izhikevich or AdEx, but where we adopt a more biomimetic form for the adaptation variable, given by the Rall's function [3]. In that model we found an unexpected analogy between the bursting state and the oscillations in a paradigmatic dynamical system of physics, the driven harmonic oscillator.
Thirdly, we showed that the problem of recursive motifs can be cast as a self-consistent problem where we provided a geometrical solution. This may open the way for progress in the understanding of spiking neurons with recursion.

In addition, from a neuromorphic electronic systems perspective, our approach has a number of innovations and advantages with respect to more conventional ones based on VLSI systems.

(1) Our synaptic circuit model introduces a novel double-step leaky integration which implements a realistic biomimetic dynamical function. It introduces a delay, which, from theoretical neuroscience, is known to be crucial to capture phenomena such as synchrony with inhibitory couplings [20].

(2) Our bottom-up methodology trades miniaturization for tunability and flexibility, hence it can be naturally made to work at biorealistic timescales. That feature is an important drawback in CMOS VLSI, whose miniature capacitors imply bioincompatible high spike rates.

(3) Our methodology provides an explicit way to avoid the mismatch problem that plagues VLSI subthreshold systems and prevents their practical adoption [15]. In contrast, our circuits are built around thyristors, thus less affected by mismatch. Moreover, using discrete components, one may achieve any desired level of device homogeneity by simply performing presorting.

(4) Our methodology is simple and of extremely low cost and availability, hence it may be immediately reproduced and widely adopted. Moreover, our SNN circuits produce robust electric signals, making it an ideal platform for the implementation of a large range of practical applications. Those systems may range from biomedical devices, such as deep brain stimulation, through reservoir computers, and neuroscience instrumentation, to central pattern generators for locomotion, and many others.

Finally, the present approach opens two exciting perspectives ahead: It may finally unlock the promise of neuromorphic electronics for artificial intelligence devices, as it provides a simple and widely accessible modular platform to build hardware devices. Also, it may provide a path to construct a general-purpose SNN simulation machine based on a sound theoretical model for research in neuroscience, beyond what can be done with software running on digital computers.

## ACKNOWLEDGMENTS

We acknowledge support from the French ANR "MoMA" project ANR-19-CE30-0020 and ANR "MemAI" project ANR-23-CE30-0040.

## APPENDIX

Here we provide details on the circuit implementations of our NS unit. We shall describe each block in detail, and give lists of materials for them along with manufacturer

<figure><img src="image_19.png" /><figcaption>FIG. 17. Details of the neurosynaptic unit circuit blocks. The neuron block is in the red box, while the blocks of two synaptic stages are in the green boxes.</figcaption></figure>

(Manuf.) in Tables I and II. We also provide the values of all components used in the different figures at the end, in Table III.

While the circuits can be easily implemented in breadboards, we also supply the printed-circuit-board (PCB) designs along with a list of components in github [62].

In Fig. 17 we reproduce the schematic circuit of Fig. 3, where we have now split the circuit into slightly different modular blocks. The convenience of this should be evident, as the two green blocks are identical. Moreover, with this partition, the system is fully modular. For instance, one may construct a network of neurons that communicate $\delta(t)$ spikes without any delay by combining the red neuron blocks. Alternatively, one may just add one green block at the output of the preneuron to implement $\exp(-t/\tau_s)$ synaptic current coupling. One may also add

<figure><img src="image_20.png" /><figcaption>FIG. 18. Schematic of the neuron plus $I_\delta$ current block. The $Q_1$, $Q_3$, and $Q_4$ are CMOS transistors from the ALD1105PBL chip. The memristor that implements the voltage-gated conductance channel of the axon hillock is implemented with a thyristor (STMicro P0118MA) with $R_m$ connected between the anode and the gate (blue box). The action-potential spike voltage $V_{ap}$ results from the fast discharge of $C_m$ on $R_a$ when the memristor is in the low-$R$ state.</figcaption></figure>
<table><caption>TABLE I. List of materials for the neuron plus $I_δ$ block.</caption><thead><tr><th>Component</th><th>Manuf.</th><th>Manuf. no.</th><th>Number</th></tr></thead><tbody><tr><td>Thyristor</td><td>STMicro</td><td>P0118MA</td><td>1</td></tr><tr><td>Capacitor</td><td>Nova</td><td>CCC-52</td><td>1</td></tr><tr><td>Resistor</td><td>Nova</td><td>CBR-11</td><td>2</td></tr><tr><td>MOSFET</td><td>ALD Inc.</td><td>ALD1105PBL</td><td>1</td></tr><tr><td>Trimmer resistor</td><td>Bourns</td><td>3296P-1-203LF</td><td>1</td></tr></tbody></table>

two green blocks, as in the schematic circuit, to implement the $α(t/τ_s)$ synaptic coupling. Another important feature, as we shall show below is that the last green block has two further options: one may implement an excitatory or an inhibitory synaptic current, both of either $exp(-t)$ or $α(t)$ type.

## 1. Neuron block

The excitatory part of the neuron circuit has been implemented in Ref. [28]; here we add the synaptic coupling. The circuit is illustrated in Fig. 18. The voltage spikes are produced on the “axon hillock” resistor $R_a$. That voltage spike needs to be transformed into a strong current impulse so one upstream neuron is able to excite, i.e., to drive spikes on a downstream neuron. Moreover, the intensity of the current needs to be modulated so as to realize the synaptic weights $W_{ij}$ between the two neurons i and j. The components to replicate the circuit are listed in Table I.

This is implemented through a simple solution, adopting a conventional current source configuration taken from

<figure><img src="image_21.png" /><figcaption>FIG. 19. The excitatory synaptic current circuit block. Here $Q_1$, $Q_2$, $Q_3$, and $Q_4$ are CMOS transistors from the ALD1105PBL chip.</figcaption></figure>

<table><caption>TABLE II. List of materials for the excitatory and inhibitory synaptic current circuit block.</caption><thead><tr><th>Component</th><th>Manuf.</th><th>Manuf. no.</th><th>Number</th></tr></thead><tbody><tr><td>MOSFET</td><td>ALD Inc.</td><td>ALD1105PBL</td><td>1</td></tr><tr><td>Capacitor</td><td>Nova</td><td>CCC-52</td><td>1</td></tr><tr><td>Resistor</td><td>Nova</td><td>CBR-11</td><td>1</td></tr></tbody></table>

the data sheet of the ALD1105PBL chip (which has two matched pairs of CMOS transistors).

The output synaptic current is $δ(t)$-like, and is implemented by a rectangular pulse of current with an intensity controlled by $R_W$. In the range of interest $I_{out} \approx 200$ µA, i.e., well within the excitatory current interval, it can be approximated by $I_{out} \approx 4V/R_W$. An important point to make is that $R_W$ can be replaced by a nonvolatile memristor in the implementation of learning neuronal networks.

The duration of the rectangular pulse is controlled by the condition $V_{spike} > \Omega \approx 0.7$ V, i.e., when the spike overcomes the synaptic threshold $\Omega$ (cf. Eq. (11) in the main text).

We can estimate some typical values to see that our neuron circuit can easily match biocompatible timescales. The duration of a spike is easily estimated, as it is directly given by the discharge timescale of $C_m$ on the axon hillock resistor $R_a$. For typical values $C_m = 1$ µF and $R_a = 1$ kΩ, we obtain $\tau_a = 1$ ms, which is a typical biological timescale.

We adopted the STMicro P0118MA thyristor, which has a trigger current $I_{GT} \approx 2$ µA. From this value and the leaky-integration time constant $\tau_m$, we can obtain an approximate expression for the firing rate for a constant input current $I_0$. The capacitor $C_m$ charges approximately

<figure><img src="image_22.png" /><figcaption>FIG. 20. The last synaptic current circuit block can be inhibitory or excitatory. For the excitatory case, $Q_1$, $Q_2$, $Q_3$, and $Q_4$ are CMOS transistors from the ALD1105PBL chip. For the inhibitory case, $Q_3$ and $Q_4$ are not used.</figcaption></figure>
<table><caption>TABLE III. Parameters for figures.</caption><thead><tr><th>Figure no.</th><th>Theoretical model parameters</th><th>Hardware model counterparts</th><th>Timescales</th></tr></thead><tbody><tr><td rowspan="3">Fig. 2</td><td>C<sub>m</sub></td><td>0.1 µF</td><td rowspan="3">τ<sub>m</sub> ≈ 68 ms</td></tr><tr><td>R<sub>m</sub></td><td>680 kΩ</td></tr><tr><td>R<sub>a</sub></td><td>2 kΩ</td></tr><tr><td rowspan="4">Fig. 4</td><td>C<sub>m</sub>, R<sub>m</sub>, R<sub>a</sub></td><td>same as Fig. 2</td><td rowspan="4">τ<sub>a</sub> ≈ 200 µs</td></tr><tr><td>R<sub>W</sub></td><td>10 kΩ</td></tr><tr><td>C<sub>s</sub></td><td>22 nF</td></tr><tr><td>R<sub>s</sub></td><td>82 kΩ</td></tr><tr><td rowspan="4">Fig. 5</td><td>C<sub>m</sub>, R<sub>m</sub>, R<sub>a</sub></td><td>same as Fig. 2</td><td rowspan="4">τ<sub>s</sub> = 1.8 ms</td></tr><tr><td>R<sub>W</sub></td><td>10 kΩ</td></tr><tr><td>C<sub>s</sub> for (a)</td><td>47 µF</td></tr><tr><td>C<sub>s</sub> for (b) and (c)</td><td>1 µF</td></tr><tr><td rowspan="3">Fig. 6</td><td>R<sub>s</sub></td><td>100 kΩ</td><td rowspan="3">(a) τ<sub>s</sub> = 4.7 s<br>(b,c) τ<sub>s</sub> = 100 ms</td></tr><tr><td>C<sub>m</sub>, R<sub>m</sub>, R<sub>a</sub></td><td>same as Fig. 2</td></tr><tr><td>R<sub>W</sub></td><td>5 kΩ</td></tr><tr><td rowspan="3">Fig. 9</td><td>C<sub>m</sub>, R<sub>m</sub>, R<sub>a</sub></td><td>same as Fig. 2</td><td rowspan="3">τ<sub>s</sub> = 470 ms</td></tr><tr><td>R<sub>W</sub>, C<sub>s</sub>, R<sub>s</sub></td><td>same as Fig. 8</td></tr><tr><td>C<sub>m</sub>, R<sub>m</sub>, R<sub>a</sub></td><td>same as Fig. 2</td></tr><tr><td rowspan="2">Fig. 10</td><td>R<sub>W</sub>, C<sub>s</sub>, R<sub>s</sub></td><td>same as Fig. 8</td><td rowspan="2">τ<sub>s</sub> = 470 ms</td></tr><tr><td>C<sub>m</sub>, R<sub>m</sub>, R<sub>a</sub></td><td>same as Fig. 2</td></tr><tr><td rowspan="6">Fig. 11</td><td>C<sub>m</sub></td><td>33 nF</td><td rowspan="6">τ<sub>m</sub> ≈ 22.4 ms</td></tr><tr><td>R<sub>m</sub></td><td>680 kΩ</td></tr><tr><td>R<sub>a</sub></td><td>2 kΩ</td></tr><tr><td>C<sub>s</sub></td><td>10 µF</td></tr><tr><td>R<sub>s</sub></td><td>100 kΩ</td></tr><tr><td>R<sub>W</sub></td><td>2.2 kΩ</td></tr><tr><td rowspan="2">Fig. 12</td><td>I<sub>0</sub></td><td>92.8 µA</td><td rowspan="2">τ<sub>a</sub> ≈ 66 µs<br>τ<sub>s</sub> = 1 s</td></tr><tr><td>C<sub>m</sub>, R<sub>m</sub>, R<sub>a</sub></td><td>same as Fig. 11</td></tr><tr><td rowspan="4">Fig. 13</td><td>C<sub>s</sub>, R<sub>s</sub>, R<sub>W</sub></td><td>same as Fig. 11</td><td rowspan="4">τ<sub>ssi</sub> = 1 s<br>τ<sub>spi</sub> = 1 ms</td></tr><tr><td>I<sub>0</sub></td><td>92.8 µA</td></tr><tr><td>C<sub>m</sub>, R<sub>m</sub>, R<sub>a</sub></td><td>same as Fig. 11</td></tr><tr><td>C<sub>s</sub>, R<sub>s</sub>, R<sub>W</sub></td><td>same as Fig. 11</td></tr><tr><td rowspan="6">Fig. 14</td><td>C<sub>m</sub>, R<sub>m</sub>, R<sub>a</sub></td><td>same as Fig. 11</td><td rowspan="6">τ<sub>ssi</sub> = 470 ms<br>τ<sub>spi</sub> = 100 ms</td></tr><tr><td>self-synapse C<sub>s</sub>, R<sub>s</sub></td><td>same as Fig. 11</td></tr><tr><td>self-synapse R<sub>W</sub></td><td>0.27 kΩ</td></tr><tr><td>projection C<sub>s</sub>, R<sub>s</sub></td><td>10 nF, 100 kΩ</td></tr><tr><td>projection R<sub>W</sub></td><td>1.5 kΩ</td></tr><tr><td>R<sub>W</sub></td><td>4.7 Ω</td></tr><tr><td rowspan="2">Fig. 15</td><td>N1 I<sub>0</sub></td><td>62.0 µA</td><td rowspan="2">τ<sub>ssi</sub> = 470 ms<br>τ<sub>spi</sub> = 100 ms</td></tr><tr><td>N2 I<sub>0</sub></td><td>46.5 µA</td></tr><tr><td rowspan="6">Fig. 15</td><td>C<sub>m</sub>, R<sub>m</sub>, R<sub>a</sub></td><td>same as Fig. 11</td><td rowspan="6">τ<sub>ssi</sub> = 470 ms<br>τ<sub>spi</sub> = 100 ms</td></tr><tr><td>self-synapse C<sub>s</sub>, R<sub>s</sub></td><td>4.7 µF, 100 kΩ</td></tr><tr><td>self-synapse R<sub>W</sub></td><td>100 Ω</td></tr><tr><td>projection C<sub>s</sub>, R<sub>s</sub></td><td>1 µF, 100 kΩ</td></tr><tr><td>projection R<sub>W</sub></td><td>0.4 Ω</td></tr><tr><td>N1 I<sub>0</sub></td><td>78.8 µA</td></tr><tr><td rowspan="2">Fig. 15</td><td>N2 I<sub>0</sub></td><td>70.6 µA</td><td rowspan="2">τ<sub>ssi</sub> = 470 ms<br>τ<sub>spi</sub> = 100 ms</td></tr><tr><td>C<sub>m</sub>, R<sub>m</sub>, R<sub>a</sub></td><td>same as Fig. 11</td></tr><tr><td rowspan="2">Fig. 15</td><td>R<sub>W</sub></td><td>2 kΩ</td><td rowspan="2">τ<sub>a</sub> ≈ 66 µs<br>τ<sub>s</sub> = 1 s</td></tr><tr><td>R<sub>a</sub></td><td>2 kΩ</td></tr></tbody></table>

Continued.
<table><caption>TABLE III. (Continued.)</caption><thead><tr><th>Figure no.</th><th>Theoretical model parameters</th><th>Hardware model counterparts</th><th>Timescales</th></tr></thead><tbody><tr><td rowspan="10">Fig. 16</td><td>C<sub>m</sub>, R<sub>m</sub>, R<sub>a</sub></td><td>same as Fig. 11</td><td></td></tr><tr><td>self-synapse C<sub>s</sub>, R<sub>s</sub></td><td>4.7 µF, 100 kΩ</td><td></td></tr><tr><td>self-synapse R<sub>W</sub></td><td>100 Ω</td><td></td></tr><tr><td>projection C<sub>s</sub>, R<sub>s</sub></td><td>1 µF, 100 kΩ</td><td></td></tr><tr><td>projection R<sub>W</sub></td><td>0.4 Ω</td><td></td></tr><tr><td>N1 I<sub>0</sub></td><td>57.3 µA</td><td></td></tr><tr><td>N2 I<sub>0</sub></td><td>52.4 µA</td><td></td></tr></tbody></table>

linearly with $V_C(t) \approx (I_0/C_m)t$, hence the gate current follows the simple expression $I_G(t) \approx V_C/R_m = (I_0/R_m C_m)t = I_0(t/\tau_m)$. The condition to fire a spike is that the gate current overcomes its switching threshold $I_G > I_{GT}$. Hence, $I_G \approx 2$ µA $= I_0$ ISI/ $\tau_m$, where ISI is the interspike interval. Thus, the firing rate results in $f(I_0) \approx (I_0/2$ µA)/$\tau_m$. For typical values of $I_0 = 60$ µA, $R_m = 1$ MΩ, and $C_m = 1$ µF, we get $f \approx 30$ Hz, which is biocompatible.

## 2. Synaptic block

The detailed circuit that implements the green blocks of the schematic model of Fig. 3 is shown in Fig. 19. The components are listed in Table II. The circuit implements the leaky integration of the input with a timescale given by $\tau_s = R_s C_s$.

This circuit equation is

$$dq/dt = -q/\tau_s + i_{in}, \quad (A1)$$

hence

$$dV_C/dt = -V_C/\tau_s + i_{in}/C_s \quad (A2)$$

and

$$\tau_s dI_R/dt = -I_R + i_{in}, \quad (A3)$$

which is exactly the form of Eqs. (2), (3), (11), and (12) of the main text. Note that, if $i_{in}(t) = I_0 \delta(t)$, then $I_R(t) \propto \exp(-t/\tau_s)$; and if $i_{in}(t) = \exp(-t/\tau_s)$, then $I_R(t) \propto \alpha(-t/\tau_s)$. This is as discussed in the main text.

The key feature that one needs to achieve with the circuit to implement the system of equations is to capture the feature of successive stage current integration. Therefore, one needs to reproduce or mirror the output current on one block as the input for the next block. As was explained in the main text, this is achieved with a standard solution, namely a current mirror. We implemented this using again the ALD1105PBL chip.

We may note that there is a price to pay for adopting this simple solution. The $R_s$ is not connected directly to the ground, as the equations above require, but to the gate

of the NMOS pair $Q_1-Q_2$ in Fig. 19. This introduces an approximation, since the approximated expression for the output current becomes $I_{R_s} \approx (V_C - 0.7V)/R_s$.

Nevertheless, as shown in the synaptic current traces in Fig. 3, the two stages are well approximated by the exponential and the alpha functional forms, which in practice demonstrates that the circuit generates the appropriate waveforms.

Finally, notice that the output current is positive or excitatory (i.e., outgoing), hence this block implements an excitatory synaptic current. The excitatory synaptic circuit can be easily modified to implement the inhibitory synaptic current case. It is sufficient to take out half of the current mirror, so the output becomes a “sink to ground” of a current of magnitude identical to its input. Hence, this would be a negative or inhibitory current; if connected to a neuron, it would decrease the charge of the membrane $C_m$. The circuit implementation is shown in Fig. 20.

With these synaptic blocks, we can realize the four model synaptic currents that we discussed in the text, as follows. For the AMPA, i.e., excitatory exponential type, adopt one block (Fig. 19). For the NMDA, i.e., excitatory alpha type, adopt two identical blocks (Fig. 19). For the GABA<sub>a</sub>, i.e., inhibitory exponential type, adopt one block (Fig. 20). For the GABA<sub>b</sub>, i.e., inhibitory alpha type, adopt two blocks, the first stage (Fig. 19) and the second stage (Fig. 20).

## 3. Model parameters

In Table III, we list all the theoretical model parameters and their counterparts in the hardware model.

[1] J. J. Hopfield, Neural networks and physical systems with emergent collective computational abilities, Proc. Natl. Acad. Sci. **79**, 2554 (1982).

[2] C. D. Schuman, S. R. Kulkarni, M. Parsa, J. P. Mitchell, and B. Kay et al., Opportunities for neuromorphic computing algorithms and applications, Nat. Comput. **2**, 10 (2022).

[3] W. Gerstner, W. Kistler, R. Naud, and L. Paninski, Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition (Cambridge University Press, Cambridge, UK, 2014).
[4] N. Brunel and M. van Rossum, Lapicque's 1907 paper: From frogs to integrate-and-fire, *Biol. Cybern.* **97**, 337 (2008).

[5] A. L. Hodgkin and A. F. Huxley, A quantitative description of membrane current and its application to conduction and excitation in nerve, *J Physiol* **117**, 500 (1952).

[6] S. R. Kulkarni, M. Parsa, J. P. Mitchell, and C. D. Schuman, Benchmarking the performance of neuromorphic and spiking neural network simulators, *Neurocomputing* **447**, 145 (2021).

[7] E. Donati and G. Indiveri, Neuromorphic bioelectronic medicine for nervous system interfaces: from neural computational primitives to medical applications, *Prog. Biomed. Eng.* **5**, 013002 (2023).

[8] C. Bartolozzi, G. Indiveri, and E. Donati, Embodied neuromorphic intelligence, *Nat. Commun.* **13**, 1024 (2022).

[9] C. Mead, Neuromorphic electronic systems, *Proc. IEEE* **78**, 1629 (1990).

[10] G. Indiveri, B. Linares-Barranco, T. J. Hamilton, A. van Schaik, R. Etienne-Cummings, T. Delbruck, S.-C. Liu, P. Dudek, P. Häfliger, S. Renaud, J. Schemmel, G. Cauwenberghs, J. Arthur, K. Hynna, F. Folowosele, S. Saïghi”, T. Serrano-Gotarredona, J. Wijekoon, Y. Wang, and K. Boahen, Neuromorphic silicon neuron circuits, *Front. Neurosci.* **5**, 73 (2011).

[11] C. S. Thakur, J. L. Molin, G. Cauwenberghs, G. Indiveri, K. Kumar, N. Qiao, J. Schemmel, R. Wang, E. Chicca, J. Olson Hasler, J.-s. Seo, S. Yu, Y. Cao, A. van Schaik, and R. Etienne-Cummings, Large-scale neuromorphic spiking array processors: A quest to mimic the brain, *Front. Neurosci.* **12**, 891 (2018).

[12] C. Bartolozzi and G. Indiveri, Synaptic dynamics in analog VLSI, *Neural Comput.* **19**, 2581 (2007).

[13] H. Wilson, Spikes, decisions, and actions: The dynamical foundations of neuroscience, (1999).

[14] S. Moradi, N. Qiao, F. Stefanini, and G. Indiveri, A scalable multicore architecture with heterogeneous memory structures for dynamic neuromorphic asynchronous processors (DYNAPs), *IEEE Trans. Biomed. Circuits Syst.* **12**, 106 (2018).

[15] S. Solinas, Dmitrii Zendrikov, and G. Indiveri, Brain-inspired methods for achieving robust computation in heterogeneous mixed-signal neuromorphic processing systems, *Neuromorph. Comput. Eng.* **3**, 034002 (2023).

[16] M. Noorman, B. K. Hulse, V. Jayaraman, S. Romani, and A. M. Hermundstad, Maintaining and updating accurate internal representations of continuous variables with a handful of neurons, *Nat. Neurosci.* **27**, 2207 (2024).

[17] E. M. Izhikevich, Neural excitability, spiking and bursting, *Int. J. Bifurcation Chaos* **10**, 1171 (2000).

[18] G. Ermentrout and D. Terman, *Mathematical Foundations of Neuroscience, Interdisciplinary Applied Mathematics* (Springer, New York, 2010).

[19] E. R. Kandel, J. D. Koester, S. H. Mack, and S. A. Siegelbaum, in *Principles of Neural Science*, 6th edn (McGraw Hill, New York, NY, 2021).

[20] C. Van Vreeswijk, L. F. Abbott, and G. Bard Ermentrout, When inhibition not excitation synchronizes neural firing, *J. Comput. Neurosci.* **1**, 313 (1994).

[21] W. Rall, Distinguishing theoretical synaptic potentials computed for different soma-dendritic distributions of synaptic input, *J. Neurophysiol.* **30**, 1138 (1967).

[22] S. Boari, G. Uribarri, A. Amador, and G. B. Mindlin, Observable for a large system of globally coupled excitable units, *Math. Comput. Appl.* **24**, 37 (2019).

[23] E. Izhikevich, *Dynamical Systems in Neuroscience* (MIT Press, Cambridge, Massachusetts, 2007), 111.

[24] S. L. Palay, C. Sotelo, A. Peters, and P. M. Orkand, The axon hillock and the initial segment, *J. Cell Biol.* **38**, 193 (1968).

[25] D. A. Wollner and W. A. Catterall, Localization of sodium channels in axon hillocks and initial segments of retinal ganglion cells, *Proc. Natl. Acad. Sci. USA* **83**, 8424 (1986).

[26] J. del Valle, J. G. Ramírez, M. J. Rozenberg, and I. K. Schuller, Challenges in materials and devices for resistive-switching-based neuromorphic computing, *J. Appl. Phys.* **124**, 211101 (2018).

[27] M. J. Rozenberg, O. Schneegans, and P. Stoliar, An ultra-compact leaky-integrate-and-fire model for building spiking neural networks, *Sci. Rep.* **9**, 11123 (2019).

[28] J. Wu, K. Wang, O. Schneegans, P. Stoliar, and M. Rozenberg, Bursting dynamics in a spiking neuron with a memristive voltage-gated channel, *Neuromorphic Comput. Eng.* **3**, 044008 (2023).

[29] H. P. Robinson and N. Kawai, Injection of digitally synthesized synaptic conductance transients to measure the integrative properties of neurons, *J. Neurosci. Methods* **49**, 157 (1993).

[30] D. Bianchi, A. Marasco, A. Limongiello, C. Marchetti, H. Marie, B. Tirozzi, and M. Migliore, On the mechanisms underlying the depolarization block in the spiking dynamics of CA1 pyramidal neurons, *J. Comput. Neurosci.* **33**, 207 (2012).

[31] M. Bikson, P. J. Hahn, J. E. Fox, and J. G. Jefferys, Depolarization block of neurons during maintenance of electrographic seizures, *J. Neurophysiol.* **90**, 2402 (2003).

[32] E. Marder and D. Bucher, Central pattern generators and the control of rhythmic movements, *Curr. Biol.* **11**, R986 (2001).

[33] E. Marder, D. Bucher, D. J. Schulz, and A. L. Taylor, Invertebrate central pattern generation moves along, *Curr. Biol.* **15**, R685 (2005).

[34] M. Ambroise, T. Levi, S. Joucla, B. Yvert, and S. Saïghi, Real-time biomimetic central pattern generators in an FPGA for hybrid experiments, *Front. Neurosci.* **7**, 215 (2013).

[35] A. J. Ijspeert, A. Crespi, D. Ryczko, and J.-M. Cabelguen, From swimming to walking with a salamander robot driven by a spinal cord model, *Science* **315**, 1416 (2007).

[36] E. Donati, F. Corradi, C. Stefanini, and G. Indiveri, A spiking implementation of the lamprey's central pattern generator in neuromorphic VLSI, 2014 IEEE Biomedical Circuits and Systems Conference (BioCAS) Proceedings, 512 (2014).

[37] H. S. Seung, D. D. Lee, B. Y. Reis, and D. W. Tank, The autapse: A simple illustration of short-term analog memory storage by tuned synaptic feedback, *J. Comput. Neurosci.* **9**, 171 (2000).
[38] J. J. Hopfield, Neurodynamics of mental exploration, Proc. Natl. Acad. Sci. **107**, 1648 (2010).

[39] K. Zhang, Representation of spatial orientation by the intrinsic dynamics of the head-direction cell ensemble: A theory, J. Neurosci. **16**, 2112 (1996).

[40] H. Seung, D. D. Lee, B. Y. Reis, and D. W. Tank, Stability of the memory of eye position in a recurrent network of conductance-based model neurons, Neuron **26**, 259 (2000).

[41] J. M. Bekkers, Synaptic transmission: Excitatory autapses find a function? Curr. Biol. **19**, R296 (2009).

[42] J. M. Bekkers, Autaptic cultures: Methods and applications, Front. Synaptic Neurosci. **12**, 18 (2020).

[43] R. Saada, N. Miller, I. Hurwitz, and A. Susswein, Autaptic muscarinic excitation underlies a plateau potential and persistent activity in a neuron of known behavioral function, Curr. Biol. **19**, 479 (2009).

[44] P. Stoliar, O. Schneegans, and M. J. Rozenberg, Implementation of a minimal recurrent spiking neural network in a solid-state device, Phys. Rev. Appl. **16**, 034030 (2021).

[45] L. Yin, R. Zheng, W. Ke, Q. He, Y. Zhang, J. Li, B. Wang, Z. Mi, Y.-s. Long, and M. J. Rasch et al., Autapses enhance bursting and coincidence detection in neocortical pyramidal cells, Nat. Commun. **9**, 4890 (2018).

[46] P. Bourdillon, L. Ren, M. Halgren, A. C. Paulk, P. Salami, I. Ulbert, D. Fabó, J.-R. King, K. M. Sjoberg, and E. N. Eskandar et al., Differential cortical layer engagement during seizure initiation and spread in humans, Nat. Commun. **15**, 5153 (2024).

[47] A. Bacci and J. R. Huguenard, Enhancement of spike-timing precision by autaptic transmission in neocortical inhibitory interneurons, Neuron **49**, 119 (2006).

[48] V. Szegedi, M. Paizs, J. Baka, P. Barzo, G. Molnar, G. Tamas, and K. Lamsa, Robust perisomatic GABAergic self-innervation inhibits basket cells in the human and mouse supragranular neocortex, *elife* **9**, e51691 (2020).

[49] C. Kittel and E. Purcell, *Berkeley Physics Course: Mechanics*, Berkeley Physics Course (McGraw-Hill, New York, New York, 1973).

[50] A. A. Hill, J. Lu, M. Masino, O. Olsen, and R. L. Calabrese, A model of a segmental oscillator in the leech heartbeat neuronal network, J. Comput. Neurosci. **10**, 281 (2001).

[51] J. E. Rubin, in *Encyclopedia of Computational Neuroscience* (Springer, New York, New York, 2022), p. 791.

[52] W. Yi, K. K. Tsang, S. K. Lam, X. Bai, J. A. Crowell, and E. A. Flores, Biological plausibility and stochasticity in scalable VO₂ active memristor neurons, Nat. Commun. **9**, 4661 (2018).

[53] C. Adda, M.-H. Lee, Y. Kalcheim, P. Salev, R. Rocco, N. M. Vargas, N. Ghazikhanian, C.-P. Li, G. Albright, and M. Rozenberg et al., Direct observation of the electrically triggered insulator-metal transition in V₃O₅ far below the transition temperature, Phys. Rev. X **12**, 011025 (2022).

[54] J. Del Valle, N. M. Vargas, R. Rocco, P. Salev, Y. Kalcheim, P. N. Lapa, C. Adda, M.-H. Lee, P. Y. Wang, and L. Fratino et al., Spatiotemporal characterization of the field-induced insulator-to-metal transition, Science **373**, 907 (2021).

[55] J. del Valle, R. Rocco, C. Domínguez, J. Fowlie, S. Gariglio, M. J. Rozenberg, and J.-M. Triscone, Dynamics of the electrically induced insulator-to-metal transition in rare-earth nickelates, Phys. Rev. B **104**, 165141 (2021).

[56] M. D. Pickett, G. Medeiros-Ribeiro, and R. S. Williams, A scalable neuristor built with Mott memristors, Nat. Mater. **12**, 114 (2013).

[57] H. Kim, M. Mahmoodi, H. Nili, and D. B. Strukov, 4K-memristor analog-grade passive crossbar circuit, Nat. Commun. **12**, 5198 (2021).

[58] C. Li, D. Belkin, Y. Li, P. Yan, M. Hu, N. Ge, H. Jiang, E. Montgomery, P. Lin, and Z. Wang et al., Efficient and self-adaptive in-situ learning in multilayer memristor neural networks, Nat. Commun. **9**, 2385 (2018).

[59] P. Stoliar, I. Akita, O. Schneegans, M. Hioki, and M. J. Rozenberg, A spiking neuron implemented in VLSI, J. Phys. Commun. **6**, 021001 (2022).

[60] N. Qiao and G. Indiveri, Analog circuits for mixed-signal neuromorphic computing architectures in 28 nm FD-SOI technology, 2017 IEEE SOI-3D-Subthreshold Microelectronics Technology Unified Conference (S3S), 1 (2017).

[61] J. Wu and M. Rozenberg, A trivial implementation of an analog spiking neuron using a memristor, for less than $1, Memristors-The Fourth Fundamental Circuit Element-Theory, Device, and Applications (2024).

[62] https://github.com/lpsssn/neurosynaptic