Characterization and modeling of spiking and bursting in experimental NbOx neuron

Hardware spiking neural networks hold the promise of realizing artificial intelligence with high energy efficiency. In this context, solid-state and scalable memristors can be used to mimic biological neuron characteristics. However, these devices show limited neuronal behaviors and have to be integrated in more complex circuits to implement the rich dynamics of biological neurons. Here we studied a NbOx memristor neuron that is capable of emulating numerous neuronal dynamics, including tonic spiking, stochastic spiking, leaky-integrate-and-fire features, spike latency, temporal integration. The device also exhibits phasic bursting, a property that has scarcely been observed and studied in solid-state nano-neurons. We show that we can reproduce and understand this particular response through simulations using non-linear dynamics. These results show that a single NbOx device is sufficient to emulate a collection of rich neuronal dynamics that paves a path forward for realizing scalable and energy-efficient neuromorphic computing paradigms.


Introduction
As the interest in Artificial Intelligence (AI) grows, spiking neural networks offer an energyefficient, hardware-compatible, and event-driven alternative to conventional artificial neural networks [1], particularly adapted for processing sensory and dynamical data.Hardware spiking neurons can be realized solely using complementary metal oxide semiconductor (CMOS) technology, but this type of implementation suffers from a lack of scalability [2].This limitation explains the growing interest in the realization of new devices that feature neuronal behavior and that can be scaled easily [3,4].However, researchers face the choice between single, scalable nanodevices that exhibit a limited range of neuronal responses and more complex neurons that offer more diverse behavior but limited scalability.Having more diverse behavior provides the potential of reproducing the brain's computational power to its full extent.Biological neurons may indeed exhibit different types of spiking responses, as well as bursting responses, where a neuron produces multiple spikes in response to an input pulse.A neuron implementing a highly simplified response will fail to provide the complexity required to emulate neurobiology.For example, the bursting response is believed to be of importance for ensuring reliable communication and synchronization between neurons [5] [6].Therefore, considerable effort has been devoted to realizing new scalable devices with diverse neuronal characteristics [5] [7] [8].
A leading idea to engineer this new type of devices is to exploit the intrinsic physics of nanoscale materials to implement neurons [9,10,11,12,13].A large number of devices have been studied for their neuronal applications [14] [15] [16]: phase change neuron [17], valence change neuron [18] [19], electrochemical metallization neuron [20], diffusive neuron [21], Mott insulator neuron [22], and spintronic neuron [23].Within these examples, metal/insulator/metal structures based on transition metal oxides such as V O x and N bO x are particularly promising candidates, as they exhibit reliable threshold switching and currentcontrolled negative differential resistance (NDR) characteristics.N bO x memristor neurons feature high endurance [24] and have been shown to be capable of leaky integrate-and-fire, all-or-nothing spiking and chaotic oscillations [25].This type of device has also been used to implement dynamic, logic, and multiplicative gain modulation [26].However, the behavior of a single device is nowhere near as complex as a real biological neuron.To obtain more sophisticated behavior, complex devices featuring multiple electrophysical processes have to be created [8], which can be challenging to model and control precisely.Alternatively, several neuronal devices can be used together in appropriately engineered circuits [7].
In this work, we fabricate and characterize memristor neurons based on a simple Pt/Nb 2 O 5 / Ti/Pt stack with current inputs and output voltage shapes that are close to the shape of a biological action potential, thanks to the effect of an inductance.These devices are straightforward to model with physics equations, and simultaneously, feature multiple computational properties such as tonic spiking, stochastic spiking, spike latency, leaky-and-fire integration (LIF), all-or-nothing firing, and phasic bursting.These neuron-like dynamics can be modelled and understood through physical equations and standard non-linear dynamics.
2 Fabrication and method N bO x memristors, comprising 5 µm × 5µm cross-point structures, were fabricated by successive film deposition and patterning.A 4-nm Ti adhesion layer and a 25-nm thick Pt layer were first deposited on a SiO 2 /Si substrate by electron-beam evaporation.These layers were subsequently patterned using optical lithography and ion-beam etching to define the bottom electrodes.A 30 nm N b 2 O 5 layer was then deposited onto the bottom electrodes using radio-frequency sputtering from a N b 2 O 5 target at room temperature in an Ar ambient.The metal-oxide-metal device was completed by adding a top electrode (10 nm Ti -25 nm Pt) deposited by electron beam evaporation.
For electrical measurements, the bottom electrode was connected to ground and the source applied to the top electrode.I-V characteristics were measured with a Keysight B1500A Semiconductor Device Analyzer after current-controlled electroforming with a positive polarity.Pulse measurements were performed using an Agilent 81160A pulse generator and a voltage-pulse to current-pulse converter (see Supplementary Materials figure 5 for the exact structure).The spiking behavior was monitored on a 2 GHz-bandwidth Keysight MSOS204A oscilloscope.All measurements were performed with a DC probe station.
Before the electroforming process, the resistance of the device was about 4 MΩ at 0.3 V. Electroforming was achieved by the application of a current ramp from 0 to 0.5 mA to the device (see Supplementary Materials figure 6).After this step, the device resistance was reduced to 93 kΩ at 0.3 V.
The simulations presented in figure 1 were computed with LTSpice using the electrical circuit presented in figure 1a based on the Newton law of cooling and the Poole-Frenkel effect (see equations ( 1) and (2) below), with a 5 ns time step.The values of all parameters used in these simulations are listed in table 1.The temperature evolution was implemented in LTSpice following guidelines described in the supporting information of [27].The simulations shown in figure 4 were executed in Python with a Runge-Kutta solver of order 5 and a timestep of 50 ps using equations ( 1) and (2) described below.

Results
The quasistatic I-V characteristics of our device are shown in figure 1a, highlighting the current-controlled S-shaped Negative Differential Resistance (NDR) response, characteristic of a voltage-controlled Threshold Switching (TS).Two characteristic values are included on the graph.The first one is the Threshold Switching point (called TS in figure 1a), where the slope of the current-controlled I-V characteristic goes from positive to negative.This point also coincides with the abrupt transition from a high-resistance state to a low-resistance state under voltage controlled transitions.The second is the hold point H, where the differential resistance becomes positive again.
The physical basis of this behavior has been under debate but is generally understood to arise from an increase in the oxide electrical conductivity due to local Joule heating.Indeed, Gibson [28] has shown that the NDR response can arise from any mechanism that gives rise to a superlinear increase in conductivity with temperature.In the case of N bO x , some authors initially attributed it to a characteristic insulator-to-metal transition (IMT) in N bO 2 [29], but it is now generally accepted that it arises from a trap-assisted transport mechanism, such as Poole-Frenkel conduction [30,31].
In the case of the Poole-Frenkel effect, a filament of oxygen vacancies connects both electrodes after electroforming.The oxygen vacancies act as potential traps for electrons.If an electric field is applied to the device, the energy profile of the conduction band in the oxide around the traps becomes asymmetric.Trapped electrons are then able to be thermally injected into the conduction band, leading to the traditional Poole-Frenkel equation for the device resistance R d as a function of the temperature T d and the voltage V d across the device: where E a is the activation energy associated with the carrier trap level, 0 the vacuum permittivity, r the relative permittivity of N bO x , q the elementary charge, and d the thickness of the oxide film.V d is the device voltage and T d is the temperature of the active device volume [30].The occurrence of electrical current through the filament results in a positive feedback, where Joule heating raises the local temperature T d , reducing the device resistance further [32,33].This phenomenon can be modeled from a lumped element model of the device, where the Newton's Cooling Law is used to describe the evolution of the temperature, where T amb is the room temperature, and C th and R th are respectively the thermal capacitor and resistor.We simulated the I-V curve of our device using these equations (see methods).
The simulation results presented with a dotted line in figure 1a show that the model reproduces the experimental data.Figure 1b presents the simple experimental setup used to measure the spiking behavior of neurons.In this circuit, R d is the device resistance described by equation 1 and C d is the intrinsic device capacitance arising from its metal/insulator/metal structure.C ext and L ext respectively account for parasitic capacitance and inductance of the measurement set-up.R out is an external resistor of 25 Ohms across which the output voltage is measured.The input of the circuit is a current, and the output is a voltage, in line with the biological configuration.Figure 1c shows an experimentally measured spike, observed by applying a constant 150 µA current input to the circuit.The shape of the output spike strongly resembles that of a biological neuron, with an initial depolarization followed by hyperpolarization: starting from a resting phase, the output voltage increases rapidly during the activation phase, and then decreases to become negative before rising again to the resting phase.
To understand this behavior, figure 1d shows simulations of the current I d flowing through the device (dots) and the simulated temperature T d of the active device volume (colored curve) during a spike, using the LTSpice model of our experiment and a current input of 180 µA (see methods).These responses are clearly correlated, with both curves exhibiting a rapid increase and a slower decrease, which can be explained as follows.The device is initially in an insulating state.When a constant current is applied, the capacitance C ext charges and the voltage across the device increases until it approaches the threshold voltage, at which point the device resistance drops, producing the increase in current and temperature evident in figure 1d.This discharges the capacitor, reducing the device voltage to the point where the memristor reverts to its subthreshold resistance.The transition to a high resistance state causes a reduction in current and temperature, ending the spike response.
The restoration part of the neuron-like voltage spike is seen in the output voltage but not in the current and temperature curves; this is due to the presence of a parasitic inductance (see figure 1b).The device intrinsic capacitance C d is small, and the current in that branch is also small.Therefore, the current going through the inductance L ext and the output resistor R out (figure 1d), is close to that going through the neuron R d .Because the voltage across the inductance opposes the variations of the current, it is first positive and then negative.The output voltage is the sum of two terms, V out = R out i out + L ext diout dt : if the inductance is large enough, the output voltage is first positive (during the activation part) and then decreases until it becomes negative (during the cooling and restoration parts).This mechanism explains the results shown in figure 1e, where the evolution of the shape of the pulse with respect to the circuit inductance is simulated.When the inductance is smaller than 100 nH, it has a negligible impact on the output voltage (showed in figure 1c); for higher inductance values, a restoration phase is observed.
Having analyzed the N bO x neuron spike shape, we now explore its computational properties.Figure 2a shows the neuron behavior when a current ramp is applied at its input.For low currents the neuron does not spike, as the NDR behavior needed for spike generation does not appear until the current reaches the TS point in figure 1a.Above this threshold current, the neuron spikes with increasing frequency until the current exceeds the hold value (H) of figure 1a, above which the NDR disappears as well as the related spiking behavior.This characteristic is reproduced in simulations in Supplementary Materials figure 7.
When the input is constant and lies between the threshold current and the hold current, the neuron spikes with a constant frequency, a behavior called tonic spiking for biological neurons, as shown in figure 2b (and reproduced in simulations in Supplementary Materials figure 8).Close to the threshold current, the behavior is stochastic, as shown in figure 2c, as can be expected from a thermally-driven process, but with a non-random occurrence of Figure 2: a: N bO x neuron output as a function of input current amplitude.A 99 µs current ramp from 0 to 0.46 mA and 1 µs fall time is applied to the device.b: Tonic spiking.The neuron receives a constant input current of 0.2 mA.c: Stochastic spiking obtained with a current of 0.109 mA.d: Spike latency.A pulse with a duration of 1 µs, a rise time and fall time of both 100 ns and an amplitude of 0.131 mA is applied to the neuron.e: Spatial integration.Comparison between two figures where a pulse of duration of 1 µs with a rise time and fall time of both 100 ns are applied to the neuron.The input current value is 0.13 mA on the left and 0.17 mA on the right.f: Temporal integration.Three pulses of of duration of 1 µs with a rise time and fall time of both 100 ns and of amplitude 0.110 mA are applied to the neuron.The frequency is 0.35 MHz on the left and 0.7 MHz on the right.spiking events, that can be described by quiet periods followed by bursts of spikes with constant frequency.Due to input current noise, the neuron output indeed fluctuates between its below-threshold behavior (no spikes) and its above-threshold behaviors (spikes with a constant frequency).This stochastic bursting behavior is reminiscent of biological neuron bursting and could be exploited for computations and learning in hardware circuits [6].
The neuron also exhibits spike latency, as evidenced in figure 2d for a 1 µs-duration pulse applied to the device.During the whole duration of the input, the output voltage does not show any significant response.However, once the pulse is back to zero, the neuron spikes.This effect can be explained naturally within the context of the above model.Indeed, when the current pulse is applied long enough for the temperature to activate the Poole-Frenkel effect, the positive feedback mechanism starts and the temperature keeps increasing even as the source stops, giving rise to spike latency.This behavior is simulated in Supplementary Materials figure 9.
Moreover, the neuron may exhibit all-or-nothing behavior.In figure 2e, two pulses with the same duration of 1 µs are applied to the neuron with different current input values: 0.13 mA for the left figure and 0.17 mA for the right one.The first pulse is not sufficient to make the neuron spike, but a slight variation of the output voltage can be observed.The second pulse is high enough to make the neuron spike, as the value of the current has been increased.In the context of a spiking neural network, this all-or-nothing behavior allows triggering a neuron only when a sufficient number of spikes (with below-threshold amplitude) arrives simultaneously at its input, thus filtering meaningful signal only, a behavior akin to spatial summation.This behavior is reproduced with simulations in figure 10.Finally, figure 2f displays a different situation where three pulses of identical duration (1µs) and peak current (0.11 mA) are applied.On the left, the input frequency of 0.35 MHz is not high enough for the neuron to spike, contrary to the right panel in which the frequency is increased to 0.7 MHz, allowing it to spike.This behavior indicates a frequency-dependent temporal summation by the neuron, reproduced with simulations in Supplementary Materials figure 11.This typical leaky-integrate-and-fire behavior is particularly adapted for spiking neural networks where frequency encodes the information.While most of the spiking features presented in figure 2 have been reported for various types of solid-state neurons [8,18,20,22], figure 3 shows that our simple N bO x neuron exhibits a behavior observed in biological neurons and scarcely investigated in memristive systems, named phasic bursting [8].In this case, for a constant input current just above the hold point (see figure 1a), the neuron starts to spike before stopping abruptly, as shown in figure 3a.This situation differs from figure 2a, where a current ramp was applied.In figure 2a, the neuron stopped spiking near the end of the input ramp, because the input current ended well above the Hold current (H point in figure 1).In figure 3a, the input is now constant and the neuron still spikes before stopping abruptly.The amplitude of the spikes appears constant, before sharply decreasing until completely disappearing.Once the neuron stops spiking, it does not start spiking again if the input does not change.Our measurements indicate that, if pulses of the right current values are applied successively, the neuron will start spiking each time before eventually stopping.However, the duration of phasic bursting is not always the same even if the input is identical.
In order to quantify the effect, a statistical study of phasic bursting as a function of input current is presented in figure 3b.A current pulse is applied to the neuron, its output is recorded on the oscilloscope, and the average frequency during the pulse duration is then computed for each point.When the phenomenon of phasic bursting occurs, spikes stop during a fraction of the total duration of the pulse, which decreases the average frequency.Despite the apparent stochastic behavior, a clear trend in the mean frequency evolution as a function of input emerges.For low currents, there is at first almost no phasic bursting, and the median frequency is almost equal to the maximum frequencies observed.Then as the input current increases, the proportion of phasic events increases and the median frequency decreases until no phasic bursting occurs.
We now present a theoretical analysis to determine the origin of the experimentallyobserved phasic bursting.We model our system with the circuit of figure 1b, neglecting the parasitic inductance and the intrinsic capacitance, that do not impact the qualitative neuron dynamics, in order to gain in simplicity and generality.The system is then simplified to two coupled first-order differential equations that link the voltage V d across the device and the temperature T d inside the active volume of the device.The first equation reads where I s is the input current and R d is the Poole-Frenkel resistance defined in equation 1.
The second equation is the Newton Law of Cooling (equation 2).Equations 2 and 3 can be solved numerically, leading to the different trajectories plotted in blue in figures 4b,c,d for the input current values I s of 0.9, 0.96702 and 1.A mA respectively.The system nullclines are also shown in dotted lines.These curves correspond to the zero values of the right-hand side of equations 3 and 2. Their intersection in the two-dimensional phase space (T d ,V d ) corresponds to points for which the derivatives of T d and V d are zero, and therefore gives the fixed point of the system for each input current.
Consistent with equation 2, the temperature nullcline does not depend on the input current I S and is therefore identical in figures 4b,c,d (orange curve).On the other hand, increasing the input current vertically shifts the voltage nullcline to the top of the phase space.The current-dependent fixed points can therefore be obtained by following the temperature nullcline.For each of these points the Poole-Frenkel resistance can be computed, and by plotting the input current I s as the function of the voltage V d (thanks to the equilibrium relation The analysis of figure 4 shows that phasic bursting is a particular situation that occurs around the hold point.Below the hold point, the fixed point is not stable, and the trajectory therefore reaches a limit cycle: this is what happens in figure 4b.At the hold point, the system undergoes a supercritical Hopf bifurcation, where the limit cycle becomes a stable equilibrium point (as seen in figure 4d).Just above the transition (figure 4c), the system reaches a stable equilibrium point, but the convergence of the trajectory is quite slow (see figure 4a).This dynamic naturally gives rise to the phasic bursting phenomenon of figure 4a, where an apparently stable train of spike unexpectedly fades out then stops.Interestingly, in the experiments, the current input range where the phasic bursting happens (∆I = 0.04 mA) is about ten times larger than in the simulations (∆I = 0.003 mA).The noise inherent to physical devices and to the input current (close to 0.018 mA in our experiments) explains the experimentally observed stochasticity of phasic bursting and expands the phasic bursting range.Indeed, even if the bias conditions of the device are set outside of the narrow range where phasic bursting is predicted in the absence of noise, fluctuations will enable the system to reach it and initiate the bifurcation, a phenomenon akin to stochastic resonance observed in biological neurons [34].Other factors can also impact the details of the phasic bursting behavior.In the model, the thermal resistance is considered constant for simplicity, but this is not true in a real device.

Conclusion
Volatile N bO x memristors are excellent neuron candidates as they are scalable, present reliable threshold switching, and are compatible with memristive synapses such as Hf O 2 Metal-Insulator-Metal structures.We have shown that the P t/N b 2 O 5 /T i/P t stack presents wellsuited I-V characteristics: a current-controlled S-shaped Negative Differential Resistance, which can be modeled by assuming Poole-Frenkel conduction.This type of devices is able to spike and the resulting shape is very close to the one of a biological neuron with initial depolarization followed by hyperpolarization due to an inductance.We demonstrated that this device presents multiple computational properties such as Leaky-Integrate-and-Fire (LIF) characteristics, all-or-nothing-firing, and phasic bursting.We also investigated the origin of phasic bursting through the analysis of the physical equations of the devices.This phenomenon comes from the bifurcation between an unstable fixed point (limit cycle) and a stable fixed point (equilibrium) driven by Poole-Frenkel dynamics.These results pave the way to easily-scalable neurons that can be easily modelled and simulated but still show a complex behavior in order to mimic biological computations.

Figure 1 :
Figure 1: a. Measured (square symbols) and simulated (dots) I-V characteristics.The V sweep and I sweep correspond respectively to the voltage-controlled and current-controlled I-V characteristics.The hold point H is indicated in green and the threshold switching point TS in red.The inset shows a sketch of the structure of the device.b. Circuit diagram of the integrated N bO x spiking neuron where C ext and L ext are respectively a parasitic capacitance and inductance.c.Measurement of a single spike of a N bO x neuron, with the four stages of an action potential indicated.d.Simulated spiking dynamics of the N bO x neuron temperature T d (colored curve) and current I d (dots) for a constant input current of 180 µA.e. Simulation of the output voltage shape with respect to the value of the circuit inductance for a constant input current of 180 µA.

Figure 3 :
Figure 3: a: Example of phasic bursting of the output voltage as a function of time.A current input of amplitude 0.47 mA is applied.The right panel zooms on the end of the phasic bursting.b: Left: Variation of the average frequency as a function of the input current.Right: Zoom on the phasic bursting regime, in order to get a statistical understanding of the phenomenon.In blue, the median frequency computed from the different average frequencies (grey dots) is plotted.

Figure 4 :
Figure 4: a: Simulations of the device current oscillations as a function of time for a current input I s of 0.96702 mA .b, c, d: simulation of the trajectory (in blue) and the nullclines (in orange for Ṫ = 0 and in green for V = 0) for different input currents I s of value 0.9, 0.96702 and 1.1 mA for each figure.The y-axis corresponds to the temperature T d in the active volume of the device while the x-axis represents the voltage of the device V d .The black arrows indicate the direction of the gradient at each point.

Figure 6 :
Figure 6: Positive current-controlled electroforming with input current going from 0 to 0.5 mA.

Figure 7 :
Figure 7: Simulation of N bO x neuron output as a function of input current amplitude.A 100 µs current ramp from 0 to 680 µA and 100 ns fall time is applied to the device.This simulation is realized in LTSpice, using the circuit shown in figure 1b and the parameters of table 1.

Figure 8 :
Figure 8: Simulation of tonic spiking.The neuron receives an constant input current of 335 µA.This simulation is realized in LTSpice, using the circuit shown in figure 1b and the parameters of table 1.

Figure 9 :
Figure 9: Spike latency.A pulse of duration of 1 µs and value 193 µA with a rise time and fall time of both 100 ns is applied to the neuron.This simulation is realized in LTSpice, using the circuit shown in figure 1b and the parameters of table 1.

Figure 10 :
Figure10: Spatial integration.Comparison between two figures where a pulse of duration of 1 µs with a rise time and fall time of both 100 ns are applied to the neuron.On the left, the value of the current is 150 µA.On the right, the input current value is 200 µA.These simulations are realized in LTSpice, using the circuit shown in figure1band the parameters of table 1.

Figure 11 :
Figure11: Temporal integration.Three pulses of of duration of 1 µs with a rise time and fall time of both 100 ns and of value 100 µA are applied to the neuron.On the left, the time period is 2.86 µs (frequency of about 0.35 MHz).On the right, the time period is 1.43 µs (frequency of about 0.7 MHz).These simulations are realized in LTSpice, using the circuit shown in figure1b and the parameters of table 1.