In order to correctly mimic biological synapses, the devices should show a controlled conductance variation. This means a modulation of the switching behavior (by means of gradual set and reset processes) to allow, in terms of ANN implementations, a regulated synaptic weight change. To do so, different voltage pulse trains can be employed. In particular, we used successive set (V_set = 0.45 V and fixed pulse widths, T_on = 1 ms, T_off = 2 ms, for a progressive set process that produces potentiation) and reset (V_reset = -0.5 V and fixed pulse widths, T_on = 1 ms, T_off = 2 ms, for a progressive reset process that leads to depression) pulse trains, as shown in Figure 2A. Multiple pulse widths and frequencies were employed in the measurements; in Figure 2, we just show the best results obtained. The voltage values employed are coherent with those found for the quasistatic I-V curves under ramped voltage stress (Figures 1C, E); in addition, they are in line with other previous works, see for instance Ismail et al. (2022). The memristive device response to successive pulse trains in terms of conductance is shown in Figure 2B for synaptic potentiation and depression.

To further demonstrate the characteristics and the reproducibility obtained with potentiation and depression stimuli, three pulse series (to allow potentiation and depression cycles) were repeated for different amplitudes (0.4 V and -0.45 V for cycle 1; 0.45 V and -0.5 V for cycle 2; 0.5 V and -0.55 V for cycle 3), while the pulse widths are fixed to 1 ms (T_on, when the pulse is active) and 2 ms (T_off, when the pulse is zero), as displayed in Figure 3A. As highlighted above, in the context of neuromorphic engineering, the pulses resemble spikes, the communication signals at the neural level. During a sequence of potentiation spikes, the memristive conductance rises. Afterward, a sequence of depression spikes leads to a conductance reduction cycle, see Figure 3B.

The device excitatory postsynaptic current (EPSC) characterizes the synaptic response to applied stimuli with different frequencies. In particular, in our EPSC study, we employed 20 Hz, 50 Hz, and 100 Hz, see Figure 4A. The stimuli consist in a train of spikes with an amplitude of 0.5 V and a time length (T_on) of 1 ms, while the time values between spikes (T_off) are 49 ms for 20 Hz, 19 ms for 50 Hz, and 9 ms for 100 Hz, see the schemes in the insets of Figure 4A. Notice that between the different spike trains, corresponding to each frequency, a 200 ms delay has been included to minimize inertial effects; after this delay time, the device operational region is assumed to cool down in what is related to thermal effects (Roldán et al., 2021). Consequently, previous signals do not affect. The EPSC increases with the pulse train frequency. This effect is depicted in Figure 4A, and it is visualized as the gain ratio of the amplitudes corresponding to the last and the first spikes in the series. The higher the stimulus frequency, the higher the EPCS gain ratio (Figure 4B). Consequently, high-frequency inputs make the synapse more active, which is beneficial for high-pass filtering in the context of spiking neural networks (Ismail et al., 2022; Li et al., 2023). At this point, it is important to highlight that in spike processing, the dynamic adaptation of the synaptic weight gives rise to many significant pattern representation and processing capabilities (He et al., 2021). In this respect, features such as EPCS are key for correctly mimicking biological synapses by means of memristors.

PPF occurs when two closely time-spaced spikes are applied to a neuron, causing the second pulse to produce a stronger response than the first. This effect is known as facilitation (Markram et al., 1997; Zucker and Regehr, 2002), and it is required for decoding temporal information in biological synapses and increasing the selectivity and information capacity of neural circuits (Zucker and Regehr, 2002). For its importance in neural processing, allowing neurons to encode data more efficiently by increasing the strength of synaptic connections between them, we have considered PPF in our analysis. We introduce two consecutive spikes (pulses generated with the semiconductor parameter analyzer, Figure 5A, in the set process operation regime of the cell), with a set delay in between, to study the corresponding synaptic response (a typical short-term synaptic plasticity effect). As explained, the first spike induces a postsynaptic response and the second induces a larger reaction. The interpulse time interval, Δt, was employed as the key variable; the shorter this interval, the higher the ratio between the average current measured for the first (I₁) and second (I₂) spikes (see Figure 5B).

Equation 1 calculates a PPF index in the usual way (Ismail et al., 2022):

Moreover, a curve can be employed to fit PPF experimental data that show an exponential dependence with the interspike time (Zucker and Regehr, 2002) (Equation 2),

where τ₁ and τ₂ are both relaxation times, and C₁ and C₂ are fitting constants. In particular, for our data (Figure 5B), the following values work correctly for the fitting: τ₁ = 0.17 ms; τ₂ = 14.42 ms, C₁ = 75.15 and C₂ = 94.13.

For our data, a simplified version of Equation 2 could work with just three parameters (τ₁, C₁ and C₂). However, the two times constant are needed if fast and slow decaying terms need to be modeled (Wang et al., 2015). See that a gradual decrease of the PPF index is obtained as the spike intervals increases. From the viewpoint of the physical mechanisms involved in the switching operation of the devices, a shorter interpulse time involves a higher temperature in the active region of the dielectric when the second spike comes in. Taking into account that the physical mechanisms behind switching are thermally activated (Dirkmann et al., 2018; Aldana et al., 2020a), the effects of the second spike in taking the set process further, and increase the device current, are more effective.

As highlighted previously, STDP is an important synaptic feature that allows the incorporation of a learning rule in spiking neural networks (Roldan et al., 2022). It can be used to implement associative learning in SNNs. Competition of spike-conducting pathways plays an essential role in establishing associations of neural connections; on the network scale, STDP potentiates the shortest neural pathways and depresses alternative longer pathways (Lobov et al., 2020). It describes the adjustment of the connection strength between neurons based on the time relation between the postsynaptic neuron and presynaptic neuron spikes in a particular synapsis (Roldan et al., 2022; Zhu et al., 2023), this mechanism is key for synaptic plasticity in biological neural circuits.

STDP characterization in memristive devices consists in the application of a delayed pair of voltage spikes to the electrodes (Roldan et al., 2022). In our experiments, the shape of the applied pulses is displayed in Figure 6A. The timing of the spikes at the top and bottom electrodes is referred to as t_pre and t_post, with the delay between them as Δt = t_post−t_pre. The bottom electrode can be left grounded to ease the measurement process and an input signal obtained subtracting the post and pre-spikes is used at the top electrode (Figures 6B–E).

In Figure 7, STDP measurements are shown. The change in device conductance (ΔG) was determined based on the starting conductance (G_INITIAL) which was obtained at the beginning of the measurement process. A good STDP behavior is obtained for different spike time widths (S_w), namely, 10, 50, and 100 μs. In order to implement the STDP as a learning rule for SNNs, Equation 3 is employed to fit the experimental data (Ismail et al., 2022; Roldan et al., 2022). A and τ parameters for potentiation and depression are employed for the experimental data fitting.

The fitting (the parameters are listed in Table 1) of the whole set of experimental data is shown in solid lines, while two other fittings to encompass the experimental dataset are depicted in dashed lines.

We have made use of the device characteristics analyzed here to build a SNN. The network architecture is shown in Figure 8E, and the operational features are given in the supplementary note 1 in the SM. The input image dataset was considered from the Modified National Institute of Standards and Technology (MNIST). The MNIST dataset is formed by 28 × 28 grayscale pixel images that consist of 70,000 handwritten digits labeled in the interval [0, 9], divided into a training set (60,000 images) and a test set (10,000 images).

We utilize the parameters shown in Table 1 coming from the fitting of Figure 7 STDP data to determine the SNN learning rule with an unsupervised learning scheme. The network input layer consists of 784 neurons, and it has been adapted to the dataset chosen in this case. Variability (as can be calculated with the constants of Table 1) was incorporated in the equation that determines the synaptic weight (Roldan et al., 2022) [traces are employed, whose value is linked to spike magnitude, the time constants in Table 1 are introduced in the differential equations corresponding to the neuron model, in our case the leaky-integrate and fire, and the A₊ and A₋ constants are employed in the equations that lead to the synaptic weight calculation (Roldan et al., 2022)]. With the new differential equation for the synaptic weight determination, including variability, we repeated the training process. Once the SNN pieces were put together, we analyzed the recognition accuracy considering a different number of epochs (Figures 8A, B) as well as a different number of neurons (Figures 8C, D). Notice that the higher the number of neurons, the better recognition accuracy for the MNIST dataset; nevertheless, the recognition accuracy improvement with the number of neurons diminish for values above 400. The inclusion of variability mostly affects the SNN accuracy with a low number of neurons; nonetheless, for 400, and mostly for 800 neurons, variability influence is low due to the SNN stochastic nature. In fact, for the higher number of neurons employed (800) and the higher number of epochs (5), there is no difference when variability is included in the calculation of the synaptic weights (see Figures 8A–D). In some experiments, higher accuracy values are obtained including variability.

The SR measurements were performed using a ramped input signal (0.28 V/s) and adding input noise with a null mean and different standard deviations (σ) (Supplementary Figure S3 in the SM). Furthermore, for the experimental SR analysis, three different statistical distributions where employed: normal or Gaussian, uniform and exponential (Heumann et al., 2016). A total of 100 I-V complete RS curves were obtained, as in Figure 1, for each standard deviation and statistical distribution. See the whole evolution of R_ON and R_OFF in the measurements in Supplementary Figure S4 (SM). There is a clear variation in the resistance evolution with rising, as expected. In particular, for the normal distribution, the variation is higher, for the exponential distribution the change in R_ON and R_OFF is found in between the results for the normal and uniform distributions. In what is connected to the set and reset voltages, notice in Supplementary Figure S5 that the added noise does not disturb much the RS operation. This result is due to the inherent stochasticity of RS operation that is resilient to added random noise. As the noise intensity rises, the difference between the set and reset voltages slightly shrinks for the three statistical distributions under study, this difference is higher for the normal distribution. As expected, the variation of the set and reset voltages increases as the noise standard deviation rises.

The cumulative distribution functions of the R_ON and R_OFF ratio are shown for different σ values and different statistical distributions in Figure 9. The CDFs shift to higher values as the noise intensity rises till approximately σ = 100 mV; at this point, the CDFs shift back. In this respect, an improvement of the device response is obtained by means of the addition of noise; in particular, at the σ value, where the resonance takes place. This behavior is clear for the normal distribution although it is not straight forward for the exponential and uniform distributions case.

The mean R_OFF/R_ON ratios vs. noise intensity was plotted in Figure 9D. A clear SR behavior is seen as it was highlighted in Mikhaylov et al. (2021) and Cirera et al. (2022). This result is in line with those shown in Rodriguez et al. (2022) although the technology employed in the study is different. In our case, SR depends on the statistical distribution function employed to generate the input noise.