Mapping Hebbian Learning Rules to Coupling Resistances for Oscillatory Neural Networks

Oscillatory Neural Network (ONN) is an emerging neuromorphic architecture with oscillators representing neurons and information encoded in oscillator's phase relations. In an ONN, oscillators are coupled with electrical elements to define the network's weights and achieve massive parallel computation. As the weights preserve the network functionality, mapping weights to coupling elements plays a crucial role in ONN performance. In this work, we investigate relaxation oscillators based on VO2 material, and we propose a methodology to map Hebbian coefficients to ONN coupling resistances, allowing a large-scale ONN design. We develop an analytical framework to map weight coefficients into coupling resistor values to analyze ONN architecture performance. We report on an ONN with 60 fully-connected oscillators that perform pattern recognition as a Hopfield Neural Network.


INTRODUCTION
Coupled oscillators have been studied for decades by scientists to describe natural phenomena (Winfree, 1967) such as the synchronization of pacemaker cells responsible for the heart beating, the synchronous behavior of insect populations, or to model neuronal activity. For instance, oscillator interactions have been shown to describe memory mechanisms and other cognitive processes in the brain (Fell and Axmacher, 2011). To characterize this variety of natural oscillations, several mathematical models (Acebrón et al., 2005;Izhikevich and Kuramoto, 2006) have been developed to explain the synchronization and phase relations in groups of coupled oscillators. Meanwhile, their massive parallel computing capability has been proved by Hoppensteadt and Izhikevich (2000), Vassilieva et al. (2011), andParihar et al. (2017) and has raised interest in designing ONN as hardware accelerators for Artificial Neural Networks (ANN) by encoding neurons' activation in the phase between oscillators. Different types of ONN have since been developed using PLLs (Hoppensteadt and Izhikevich, 2000) or oscillators in CMOS technology (Maffezzoni et al., 2015a(Maffezzoni et al., , 2016Jackson et al., 2018), demonstrating pattern recognition or resolution of some optimization problems like the Traveling-Salesman-Problem (Endo and Takeyama, 1992). However, to design a competitive ONN at a large scale, a design framework is needed to establish a formalism on how to perform computations with ONN and also compare its energy efficiency with ANNs running on digital processors.
Researchers have developed oscillators by using non-linear devices such as spin-torqueoscillators (Raychowdhury et al., 2019) or with materials presenting a hysteresis resistive state to induce electrical oscillations when properly biased (Sharma et al., 2015;Wang et al., 2017). A compact device that transitions from metallic to insulating state (MIT) can be manufactured with vanadium dioxide (VO 2 ) (Corti et al., 2018), and has recently become an interesting candidate to design energy-efficient relaxation oscillators. Moreover, coupled-VO 2 -based oscillators have been experimentally validated for various applications such as image saliency detection (Shukla et al., 2014), graph coloring (Parihar et al., 2017), filters in Convolutional Neural Networks  and implementing Hopfield Neural Networks (HNN) for pattern recognition (Corti et al., 2020).
In an HNN defined by Hopfield (1982), every neuron is connected to all the others, and synaptic weights are computed with the Hebbian rule (Hoppensteadt and Izhikevich, 2000). However, setting the right coupling element between oscillators in ONNs remains a challenge . Given N fully-coupled VO 2 -oscillators, it is yet unknown how to transform the coefficients obtained analytically via the Hebbian learning rule to coupling resistor values among oscillators. Further, how can one interpret the coefficient signs W ij such as positive or negative values?
For weakly coupled oscillators with sinusoidal waveform, one can use the models that exist in literature (Izhikevich and Kuramoto, 2006;Maffezzoni et al., 2016) for synaptic design. However, it is more difficult for non-linear relaxation oscillators as there is no direct mapping between models and hardware. Two VO 2 -oscillators coupled by a capacitance or a resistance have been studied (Maffezzoni et al., 2015b;Parihar et al., 2015), but to the best of our knowledge, there is not yet a formalism to map weights to coupling elements in a larger network. This formalism is a crucial step to allow largescale ONN design exploration. A greedy approach would be to tune the coupling elements corresponding to the most negative and most positive weights and linearly interpolate all the other weight values. However, this would be impractical. It would require repeated simulations and re-tuning coupling resistances when changing any oscillator parameter; hence, it is not suitable for large-scale ONN design. Parihar et al. (2015) proposed to use capacitors or resistors to implement a negative or a positive weight, respectively. However, it would imply using twice as many components to emulate a complete signed synaptic range.
In this work, we propose a mathematical framework to map both negative and positive Hebbian coefficient values to ONN coupling resistances, as illustrated on Figure 1. We first present a single VO 2 -oscillator followed by the dynamics of two coupled oscillators. Then, we show that adding switches between oscillators and coupling elements enhances the ONN dynamics control. Based on this simple architecture, we present the ONN computation style and how coupling resistances set the ONN memory expressed in different phase states. Next, by merging oscillators' dynamics with HNN formalism, we introduce the ONN mapping function that maps Hebbian coefficient values to coupling resistance values. Finally, we report on the architecture and mapping results by simulating 60 coupled VO 2 -oscillators for pattern recognition.

General Properties
Unlike most ANNs where signals of interest are amplitudes, ONN consists of N VO 2 -oscillators coupled by resistances where the final result is given by N-1 phase relations to a reference oscillator (Corti et al., 2018). Despite this difference, there are several aspects common to any ANN that motivate our ONN study: • The network is composed of neurons (oscillators) having input and output nodes. • The activation function between the oscillator input and output is non-linear. For more than two oscillators, the activation function is unknown but bounded as the output phase difference of neuron i φ out i is in the range [0; 180 • ]. • ONN has a memory (coupling resistances) and can therefore be trained to achieve specific functionality.
We use these properties to implement an HNN (Hopfield, 1982) by encoding its binary outputs in the ONN phase as shown in Figure 1A. Conveniently, we can represent each oscillator's output as a black pixel if φ out i = 180 • or as a white pixel if φ out i = 0 • .

VO 2 Device
Vanadium dioxide is a material which presents temperaturedriven phase change transitions. At room temperature, it remains semiconductor (insulating state) with a monoclinic crystal structure and transitions to a rutile metallic state when the temperature reaches a threshold . VO 2 presents an hysteresis behavior as it needs to reach a lower temperature threshold to transition back to the insulating state (Figure 2A).
When a VO 2 device is in series with a biasing load (Figure 2B), the voltage drop V across the device induces Joule heating which can trigger the insulating to metallic transition (IMT). To obtain the VO 2 IV characteristic, we sweep the supply voltage which triggers IMT and MIT when V ≥ V H and V ≤ FIGURE 2 | (A) VO 2 resistance vs. temperature. R VO2 axis is in logarithmic scale. (B) VO 2 oscillator circuit (C) VO 2 I-V curve showing the device hysteresis behavior. When its voltage reaches a threshold V H , it transitions from insulating to the metallic state. To transition from a metallic to insulating state, V must be lower than the threshold V L . (D) If the device is biased in its Negative Differential Resistance region (NDR), its state alternates between the metallic and insulating state, producing electrical oscillations.
V L , respectively. The VO 2 hysteresis behavior appears in the IV characteristic with a typical Negative Differential Region (NDR) suitable to bias the device and produce oscillations ( Figure 2C).
In this work, we use the VO 2 compact model from Maffezzoni et al. (2015b) which reproduces the VO 2 hysteresis along with continuous and abrupt transitions between the two states. The VO 2 hysteresis behavior is conceptually emulated by an amplifier with positive feedback that charges or discharges a RC circuit when the VO 2 is in metallic or insulating state, respectively. The voltage V c across the capacitor commands the VO 2 conductance G VO 2 as: Where R ins and R met are the VO 2 resistances in insulating and metallic state, respectively. The dynamics of the VO 2 conductance is given by the RC circuit where τ 0 is its time constant modeling the transition time of the VO 2 : V 0 (t) is the output of the positive feedback amplifier (gain α) and is expressed as:

VO 2 Oscillator Circuit and Dynamics
We bias the VO 2 device with a series resistor R S to obtain a compact relaxation oscillator. To produce oscillations, the load line I L set by V DD and R S must intercept the VO 2 I-V curve in its NDR to obtain an unstable fixed point ( Figure 2C). The VO 2 T osc 21.6 µs device state hence alternates between the metallic and insulating state. When the VO 2 device is in the insulating state, the parallel capacitance C P at the output node discharges through R S until the VO 2 voltage reaches V H and transitions to the metallic state. Then, C P charges through the VO 2 device until its voltage reaches V L and a new cycle begins ( Figure 2D). We use circuit parameters depicted in Table 1. In this case, the load resistance R S is 20 times larger than the VO 2 metallic resistance R met ; thus, the capacitance charge time is much faster than its discharge ( Figure 2D).
Oscillator's dynamics can be described by Kirchoff 's law as: Note that despite the first order differential equation, oscillations can occur as Equations (1-3) describe the hysteresis behavior of G VO 2 (t). To solve the oscillator dynamics, we start from an initial insulating VO 2 state and solve numerically on Matlab the system of Equations (1-4) by using Euler's method and Newton-Raphson's algorithm for non-linear Equation (3). Figure 2D shows an example where V out (t = 0) = 0 V, G VO 2 (t = 0) = 1/R ins , and V DD (t = 0) = 2.5 V.

Initialization of Two Coupled Oscillators
Two coupled oscillators represent the smallest ONN and serve as the building block for large-scale ONN. To provide input to the ONN, we delay the second oscillator V DD starting time with respect to the first oscillator (reference oscillator) to set an initial phase relation between them. Assuming oscillators have the same period T osc , we can translate the input delay t init as an initial phase relation as: However, if the two oscillators are always connected, they might have different oscillation periods during initialization. Therefore, their initial phase relation cannot be represented as a proportion of T osc (5). For example, as shown in Figure 3A, the second oscillator starts t init = 0.5 T osc after the first one to set an FIGURE 3 | (A) Two oscillators are coupled by a resistance R C =10 k without coupling switches in between. V DD2 is turned-on 0.5T osc after V DD1 to set an initial phase of 180 • . However, for t <0.5T osc the first oscillator period is decreased due to the shunt R C at its output node, and we cannot control the initial phase difference. The two oscillators are in-phase after convergence. (B) Two oscillators coupled with resistance R C =10 k and coupling switches between isolate the oscillators during initialization. This time, the 0.5T osc input delay sets desired 180 • initial phase state. Here, the switches are closed at t on = 0.5T osc + t c such that V out2 (t on ) = V + . The two oscillators converge to an 180 • phase state relation.
initial phase relation φ init = π. For t < t init , the second oscillator is off, and its output node is floating. Therefore, during this time, the equivalent load resistance of the first oscillator is R S //R C , which induces a shorter period of oscillation T ′ osc < T osc and hence no control on the initial phase.
We introduce switches between each oscillator and coupling elements as in Figure 3B to tackle this lack of control. We let each oscillator switch freely with a known oscillation period T osc before coupling them at t on . Their dynamics can be expressed by Equation (4) and the initial conditions will be known at t on . ONN initialization is improved at the cost of one additional switch per oscillator, which can be achieved with a transfer gate.

Dynamics of Two Coupled Oscillators
To predict the output phases and demonstrate ONN ability to store information, we express the dynamics of the two coupled oscillators using Kirchhoff 's laws: Currents are I c1 = −I c2 representing the coupling element's current flow. As for the single oscillator case, we numerically solve (Equation 6) along with VO 2 (Equations 1-3). Figure 4 shows a simulation where V DD2 is turned on 0.1T osc after V DD1 which initializes a light-gray pixel for oscillator 2 input image. For a small coupling resistance, R C =10 k , ONN converges to a stable state with both oscillators in-phase (0 • ,0 • ). Whereas, for R C =100 k , ONN converges to out-of-phase (0 • ,180 • ). In the next subsection, we study the role of R C on ONN memory and investigate how to retrieve a stored pattern by applying an input delay t init . This formulation is the core of our proposed mapping function to translate Hebbian coefficients to ONN coupling resistances.

Memory of Two Coupled Oscillators
We solve numerically (6) and extract the output phase relation between oscillators. Figure 5A shows the simulation results. As already observed by Corti et al. (2018), a large coupling resistance R C > 40 k induces oscillators in out-of-phase relation (0 • , 180 • ) for any input delay, whereas a small coupling resistance R C < 10k aligns oscillators in-phase (0 • , 0 • ) for any input delay.
In contrast, we examine the region between these two ranges, highlighted in the center of Figure 5A. We observe that for 10k ≤ R C ≤ 40k both states co-exist and oscillators store two patterns (0 • , 0 • ) and (0 • , 180 • ) that can be retrieved by adjusting the input delay. The line transition function between in-phase and out-of-phase regions represents our analytical function for ONN memory with respect to coupling resistance and initialization. It is defined as: Frontiers in Neuroscience | www.frontiersin.org With t transit the initial delay such that: To confirm the existence of ζ (R C ), we emulate VO 2 oscillators with off-the-shelf components on a Printed Circuit Board (PCB) and we reproduce the experiment of two coupled oscillators. The relaxation oscillator circuit consists of an inverting Schmitt trigger (Schmitt, 1938) Operational Amplifier (OPA) that implements the VO 2 hysteresis behavior ( Figure 5B). The OPA saturates to +V sat and −V sat while the 1.8 nF output capacitor charges and discharges, respectively. Figure 5C shows the voltage across the output capacitor for a decoupled oscillator. Similarly to a VO 2 oscillator, the OPA transitions to another state when the voltage across the output capacitor reaches V + or V − . The 5.6 k resistor implements the metallic VO 2 resistance, whereas the 100 k resistor corresponds to the load R S . For a fixed oscillating period of T osc =200 µs, we vary R C and we measure t transit values that define the experimental transition function ζ (R C ) ( Figure 5D). There is a good match between experimental ζ (R C ) data points and the analytical transition function derived in next subsection. Such formulation ζ (R C ) is of interest as it represents a closed-form representation of ONN memory instead of repeating numerical simulations for different oscillator parameters.

Phase Transition Function for Two Coupled Oscillators
The phase transition function has already been observed (Nez et al., 2021) but to the best of our knowledge, no closed-form expression has ever been reported. To obtain the transition function, we solve node voltage equations analytically for two coupled oscillators during initialization (see Supplementary Material). We derive oscillator outputs as: V 0 out1 and V 0 out2 are the initial voltages when oscillators are coupled and τ ′ is defined in Supplementary Material Equation (S18). Equation (9) describes both oscillator output voltages attracted via the coupling resistance R C . If the coupling is strong enough (small R C ), both oscillators are rapidly pulled together with a speed determined by τ ′ . If V < ǫ (ǫ defined in Supplementary Material Equation S22) before reaching the VO 2 threshold V − , then both oscillators will transition to low resistive states, and the exponential term in Equation (9) will keep the two voltages locked. This concept is illustrated in Figure 6A when both oscillators are in-phase. However, if V out1 reaches V − before V out2 such that V > ǫ as in Figure 6B, the first oscillator transitions to a low resistance state (metallic state) while the other oscillator is still in high resistance state (insulating state). The two oscillators are then in opposite states, and this leads to out-of-phase relation.
Thus, to obtain the transition function ζ (Equation 7), we study the case when both V = ǫ and V out1 = V − conditions are fulfilled (see Supplementary Material for details). By combining (Equations S17, S24, and S25 in Supplementary Material) we derive coupling resistance as:  (10), and obtain a relation between R C and t transit . Note that ǫ is a function of R C , thus we cannot solve (10) analytically. Instead, we numerically solve (Equation 10) using Newton-Raphson's algorithm for t transit values. We finally obtain R C values that describe the inverse of the transition function as: Transition function ζ is plotted as the curve line (red line) in Figure 5A, and there is an excellent fit between our analytical model, simulations (transition region between dark and light green in (R C , t init ) plan) and the transition curve obtained experimentally with off-the-shelf relaxation oscillators ( Figure 5D). In addition, we can now extract the coupling resistor R 0 that corresponds to a neutral synaptic connection W = 0. As by definition, both output phase states can equally occur for W = 0, we extract R 0 as: FIGURE 5 | (A) Plot showing the phase relation between two oscillators for every set of parameters (R C , t init ). As expected, small coupling resistances tend to pull the oscillator phase together, whereas large coupling resistances push the phase away. The green region shows the coupling resistance range 10k < R C < 40k in which two patterns (0 • and 0 • ) and (0 • and 180 • ) are memorized and can be retrieved by adjusting the input delay. The red curve is our analytical model describing the transition between the two phase states in the plan (R C , t init ), and plays a major role in the ONN ability to memorize patterns. (B) Experimental set-up of two coupled relaxation oscillators based on MCP6001 OPAs. We delay VG2 with respect to VG1 by t init to set the initial phase, and we close SW after initialization. Finally, based on the transition function, we predict the final phase relation as: Analogous to ANNs, Equation (13) can be thought of as oscillator's activation function. Because, it provides the oscillator's output phase based on its input phase (set by t init ; Equation 5) and the weight implemented by R C .

Impact of VO 2 Parameters Variations on the Phase Transition Function
Fabricating reliable VO 2 devices is challenging (Corti et al., 2019) and ONN experiments with VO 2 are currently limited to few devices because of device variability (Shukla et al., 2016).
Here, we study the impact of VO 2 variability on the ability to phase-lock and on the synaptic range. The transition function ζ (R C ) defines the boundary between two phase regions, and allows direct identification of the neutral coupling resistor R 0 corresponding to the weight W = 0 (Equation 12). Thus, we use ζ and R 0 as metrics to assess the impact of VO 2 parameters' variations. We apply relative variations on VO 2 parameters one at a time from -20% up to +20%, as shown in Figure 7A. Note that we vary V H from -4% up to +4% only, as for larger positive variations oscillations do not occur (load line crosses the insulating branch and forms a fixed point). For all cases, we discretize the whole input space (R C , t init ) and perform multiples transient simulations to extract the new phase regions. Then, we numerically solve (Equation 10) with the new sets of parameters and we verify that the transition function ζ matches the phases boundary obtained via transient simulations, as in Figure 5A. Figure 7A shows the set of phase transition curves obtained when varying R met , R ins , V L and V H . Note that our current formalism assumes matched oscillators and hence, variations are FIGURE 6 | Two identical oscillators coupled by R C = 12k . (A) The second oscillator is turned on at 0.2 T osc after the first one. By zooming on the waveform when the first oscillator reaches V − , we observe a voltage difference V < ǫ. Therefore, the oscillators converge to an in-phase state. (B) The second oscillator is turned-on 0.3 T osc after the first one. In this case, we observe a voltage difference V > ǫ and the oscillators converge to an out-of-phase state.
applied to both coupled oscillators. For all curves, the maximum t init value corresponds to an input delay of T osc /2 and shows the oscillation period variation (highlighted in green in Figure 7A). Finally, we extract R 0 for each configuration ( Figure 7B). We observe that variations on the IMT point (defined by R ins and V H ) induce the largest T osc and R 0 variations. With our biasing set by R S and V DD (Table 1), the most sensitive VO 2 parameter is V H as +4 and -4% V H variations induces +40 and -20% R 0 variations, respectively. As the dynamic of the voltage V across the VO 2 device is given by C P dV/dt = I L − I, we believe this sensitivity is mainly due to the load line that passes very close to the IMT point on the VO 2 IV characteristic ( Figure 7C). In this case near IMT, I L (V H ) − I(V H ) is small and the voltage "slows down" and is very sensitive to any IMT variation. When applying -4% up to +4% V H variations, the oscillating period T osc almost doubles (same remark with -20 and +20% R ins variations). Ideally, we would then place the load line at equal distances between MIT and IMT points (I(V L ) − I L (V L ) ≈ I L (V H ) − I(V H )) to homogenize the impact of VO 2 variations. However, we show in the next subsection that such biasing would prevent binary phase locking and that resistively coupled oscillators need a very asymmetric waveform to phaselock to 180 • .

Impact of Oscillators' Waveform Shape on ONN Phase-Locking
Oscillators' circuit parameters listed in Table 1 influence the oscillating frequency, amplitude and waveform shape. The oscillating waveform shape has a major influence on ONN phase-locking capability and has been studied for PLL-based ONNs by Hoppensteadt and Izhikevich (2000). Here, we study the impact of the oscillating waveform shape on the capability for pairs of oscillators to lock to the 180 • phase state. We characterize the oscillating waveform shape with the ratio τ d /τ c , where τ d and τ c are the discharging and charging time constant, respectively (defined in Supplementary Material Equations S6, S7). Our transition function ζ links ONN phaselocking properties to the metric τ d /τ c , as ζ only depends on oscillators' internal parameters.
We reproduce the previous simulation with two coupled oscillators to extract the output phase regions for different load resistances R S that set τ d /τ c (Figure 8A). Note that we also could have varied VO 2 parameters such as R met , but instead we consider the same device. For τ d /τ c = 3.7 (R S = 3 k ), we observe that the two oscillators cannot lock to φ out = 180 • for small t init values. In other words, the phase state φ out = 180 • stored by a large R C cannot be fully recovered. This can be an issue for some FIGURE 7 | (A) Phase transition curves ζ (R C ) when varying VO 2 parameters R met , R ins , V L from -20 to +20%, and V H from -4 to +4% for circuit parameters listed in Table 1. Left hand side of the transition curve corresponds to inputs ( t init , R C ) inducing φ out = 0 • whereas right hand side corresponds to φ out = 180 • phase region. V H and R ins variations correspond to IMT point variations and have the most detrimental impact on the transition function variations. The oscillating period almost doubles due to IMT variations. (B) Variations of the neutral synaptic resistance R 0 with respect to VO 2 parameters' variations. R 0 is very sensitive to V H as -4 and +4% V H variations induce -20 and +40% R 0 variations, respectively. (C) The ONN sensitivity to IMT point is mainly due to the load line I L = (V DD − V)/R S placed close to IMT point. Any IMT variation greatly impacts the oscillators' dynamics defined by C P dV/dt = I L − I. pairs of oscillators that need an out-of-phase relationship for any input delay.
If τ d /τ c = 59 (R S = 20 k ), the charging time is much smaller than the discharging time and the oscillating waveform becomes very asymmetrical. Interestingly, this configuration enlarges the 180 • phase region and φ out = 180 • is reachable for any t init value for large R C . Our analytical model ζ (R C ) predicts the correct boundary between the two phase regions (red plain lines in Figure 8A).
We study a simple case where 4 VO 2 -oscillators are coupled by resistances to store a single pattern ( Figure 8B). Based on transition functions obtained for 2 coupled oscillators, we compute coupling resistances R +1 and R −1 that correspond to synaptic coefficients +1 and -1, respectively. We set R +1 and R −1 around R 0 as R +1 = ζ −1 (T osc /4 + T osc /8) and R −1 = ζ −1 (T osc /4 − T osc /8), respectively. Then, we scale coupling resistances as 3xR +1 and 3xR −1 as every oscillator is connected to 3 others ( Figure 8B). We notice that ONN with τ d /τ c = 59 retrieves the correct stored pattern whereas ONN with τ d /τ c = 3.7 produces a wrong output (Figure 8C). In the latter case, we observe that all oscillators converge to an in-phase relationship. We believe this wrong behavior is mainly due to the small τ d /τ c value for which it is less likely ONN converges to φ out = 180 • , as described by the transition function ζ . Figure 9 shows results of the same experiment for τ d /τ c varied from 1.8 up to 59 (obtained for 2 k ≤ R S ≤ 20 k ). We observe that τ d /τ c > 20 is required to retrieve the correct pattern.
Interestingly for τ d /τ c ≤ 20, there are cases where the fourth oscillator locks to a phase state around 270 • . 270 • phase value is also obtained in the phase plot between two coupled oscillators, such as on the left-hand side of Figure 8A. This phenomenon would allow more than two phase values but is not captured by our current formalism. In contrast, we set τ d /τ c to high values (59 in this work) to ensure binary 0 • and 180 • phase locking.
To implement a large-scale neural network such as HNN with an ONN, we need a systematic approach to map the weights to ONN coupling resistances. In the next section, we exploit some HNN features and our knowledge of two coupled oscillators to propose a mapping function.

Applying HNN Formalism to ONN
We exploit HNN formalism to build an analogous representation in ONN. For equivalence, we treat HNN neurons similar to ONN oscillators. Such as, we consider a neuron i with two possible states S i that can be thought of as equivalent to ONN oscillators with 0 • or 180 • phase relations as: In HNN, each neuron output state is dynamically determined by a sigmoid activation function g(x) = tanh(βx) + 1 /2 (with β a positive parameter) (Gerstner et al., 2014) and shown in FIGURE 8 | (A) Phase plots showing φ out with respect to t init and coupling resistance R C between two oscillators. For τ d /τ c = 3.7, 180 • phase state is not reachable for low t init values. In contrast for τ d /τ c = 59, 180 • phase-locking can occur for any t init value for large R C . The red line is our analytical model ζ (R C ) and captures well the boundary between phase regions. (B) Four coupled oscillators store a pattern composed of 2 white and 2 black pixels. Positive and negative weights are mapped to 3xR +1 and 3xR −1 , respectively. (C) ONN inference for τ d /τ c = 3.7 and τ d /τ c = 59. The first configuration leads to a wrong in-phase relationship for all oscillators. In the latter case, ONN retrieves the correct stored pattern. Figure 10. For a neuron i, g gives the probability to reach one of the two states at t + t for a given input weighted sum h i (t) as with In ONNs, Equation (15) would represent the probability of oscillator i to be in-phase with the reference at time-step t.
For two oscillators case, the weighted input sum of the second oscillator is given by: Then, the probability of the second oscillator to be in-phase with the reference can be derived by Equations (15) and (17), as:

Mapping Function
Here, we apply the above definitions to derive a mapping function using the HNN formalism, as: where, normalized weights are −1 ≤ W ij ≤ 1 and N is ONN size. Before scaling to N oscillators, we derive a mapping function µ 2 for two coupled oscillators. The unifying step between HNN and ONN is the recasting of the phase transition curve, ζ as the probability P inphase for a given coupling resistance R C . In ONNs, the input delay t init can be considered as a uniform random variable taking values between 0 and T osc /2 and the transition function ζ would give the probability P inphase (for example for R C > 10k ): and by Equations (18) and (20), we finally obtain This mapping function is represented in Figure 11 for three different values of the sigmoid parameter, β. We see that this (C) Output phase with respect to τ d /τ c . For τ d /τ c < 7, all oscillators converge to a wrong in-phase relationship. For 7 ≤ τ d /τ c < 20, the fourth oscillator locks to a 270 • phase state. A very asymmetrical waveform such that τ d /τ c ≥ 20 leads to a correct binary phase-locking.
FIGURE 10 | Model of artificial neuron used to construct our mapping function µ N . The neuron's output state S i (t) is either +1 or -1 and is dynamically updated at each time-step t according to the sigmoid activation function g.
Here, g gives the probability to have one of the two states at t + t for a given input weighted sum h i (t).
parameter sets the range of R C and could be adapted for different ONN sizes. Interestingly, we notice that | W 21 / R C | is quite large for a positive weight, whereas it is much smaller for a negative one. For example, we see in Figure 11B that the function ζ −1 is a logarithmic function; thus, any small variation in R C around 10k is likely to change the final phase state outcome. Whereas, for large R C , the two oscillators are almost always out-of-phase. This asymmetry in ζ −1 comes from the oscillator waveform type, as ζ −1 is derived from Equation (10), which is specific for relaxation oscillator waveform type. Hence, we expect some change for other types of waveforms, such as linear sawtooth, but the formulation of mapping (21) is general enough to be applied to any relaxation oscillators. For large-scale ONN with N oscillators, we scale µ 2 (21) by a factor N − 1 to ensure the conservation of the current flow in coupling resistances. We finally obtain: In next section, we demonstrate the effectiveness of the proposed mapping function (22) to design a 60-ONN architecture for pattern recognition as in HNN.

ONN Training and Mapping
In the previous section, we presented the memory capability of two coupled oscillators. Here, we apply the analytical formulas to a larger ONN size. We develop a design flow as shown in Figure 12A for pattern recognition with ONNs where we have implemented the proposed mapping function. We first compute the weights associated with the M stored patterns using the Hebbian Rule (Hoppensteadt and Izhikevich, 2000), as: We store M = 6 images representing digits "0", "1" to "5" as shown in Figure 12B. Next, we use our mapping function to compute the coupling resistances associated with the Hebbian coefficients. The mapping is represented in Figure 12C for different values of parameter β which sets the slope of µ N (W ij ). Increasing β induces a larger coupling resistance range. Because the Hebbian rule normalizes the weights by the network's size N (23), we scale β with N to keep a relative standard deviation of R ij approximately constant when increasing the ONN size. By simulations, we found that the best accuracy is obtained for β = N/32, and we report the results in subsection 3.1.3.

ONN Inference
For every input, we set up the 60 ONN with black and white pixels encoded by -1 and +1, respectively. For pixels with black input, the corresponding oscillator is initialized with a delay t init = T osc /2 to set an initial out-of-phase relation (5) whereas, for a white pixel, no delay is introduced. A noisy gray pixel corresponds to an input delay between 0 and T osc /2. After few oscillations, the ONN settles, retrieves the noiseless pattern and phase relations are measured. An example of ONN voltage dynamics is presented in Figure 13A, showing the initialization and the settling time before the ONN stabilizes. Figures 13B,C show two examples of input images where 15 pixels have been randomly altered by a uniform distribution taking values between -1 and +1. When the number of noisy input pixels is too large such as in Figure 13D (20 noisy pixels), ONN converges toward a wrong spurious state that is different from the stored patterns. The results are in accordance with original observations from Hopfield (1982), proving that our mapping can implement HNN with ONN. Inverse of the transition function ζ −1 determines the coupling resistance R C for a given probability P inphase . (C) The mapping function µ 2 is obtained by the composite function ζ −1 (g).

FIGURE 12 | (A)
Illustration of the ONN design flow for the associative memory application. Patterns to store can be represented as black and white images from which we compute weights with the Hebbian rule during the training process. Then, we use the mapping function µ N to get the coupling resistances, allowing a systematic ONN design. (B) Stored patterns. (C) Coupling resistances as a function of Hebbian weights, computed with the mapping function µ N for different values of parameter β.

ONN Recognition Accuracy
Here, we perform simulations to compute the pattern recognition accuracy of the 60-ONN. We randomly apply noise to training patterns to generate a test set. It consists of 20 different subsets S k , k ǫ {1, 2, .., 20} in which 60 different test patterns have k randomly located fuzzy input pixels. We vary the mapping function parameter β to assess its influence on ONN accuracy. We notice from Figure 14C that ONN achieves the best accuracy for an optimum value β = N/32 = 1.875. In this case, ONN recognizes more than 80% of test images with up to 20% of FIGURE 13 | (A) Noisy input image "2" with 15 random altered pixels and voltage waveforms of 60 oscillators. ONN is initialized during T osc /2 with the noisy input image. After few oscillation cycles, ONN settles, and phases are measured. ONN retrieves the correct output image that corresponds to the stored pattern "2." (B) The input image "5" has been altered at 15 random pixel locations by a uniform distribution. ONN retrieves the corresponding stored pattern. (C) Similarly, 15 random pixels of the input image "2" are altered, and ONN retrieves the correct corresponding pattern. (D) In this example, 20 noisy input pixels are introduced to digit "1," and ONN converges toward a spurious state.
noise. As seen in Figure 12C, the slope parameter β sets the coupling resistance range, which in turn affects ONN accuracy. For instance, we observe that the set of coupling resistances obtained for β = 2 is similar to the case β = 2.5, but the accuracy is much lower in the latter case. In the next section, we quantify the coupling resistance accuracy that is required for synaptic design.

ONN Coupling Resistance Range
We study the impact of R C 's relative variations and R C 's mean value. R min C is the minimum resistance common to all coupling resistances, and R is the additional series resistance to distinguish between weights ( Figure 14A). Using the Hebbian rule, weights are located near "0" coefficient as in Figure 12C and our mapping function can be fitted linearly (dashed lines). Therefore, every coupling resistances can be approximated by R C ≈ R min C + n R min with n ǫ {0, 1, 2, .., M}. Using this linear approximation, we can verify ONN accuracy is similar to the nominal case of mapping function with β = N/32, as shown in Figure 14C with the magenta dashed line.
As observed in previous sections, ONN accuracy is quite sensitive to the coupling resistances. We obtain R max ≈ 15% R min C for β = 3, and R max ≈ 5% R min C for β = 1 (Figure 14B). For these two cases as shown in Figure 14C, ONN shows poor accuracy. It is rather for β = N/32 = 1.875 with R max ≈ 10% R min C that ONN accuracy is above 90%.
To achieve the best ONN accuracy, a very good resistor matching is required, as we need a precision of R min = 1.7% R min C between two consecutive weights. To study the influence of the R C 's mean value only, we apply the same variation to all coupling resistances for β = N/32 and for 10 fuzzy input pixels ( Figure 14D). We notice that the mean value of coupling resistances can vary from -10% up to +5% from its nominal value to achieve a similar accuracy.

DISCUSSION
Oscillatory neural networks are triggering great interest for parallel processing applications, but a remaining challenge is how to compute with ONNs. To do so, we build analogies with ANN to determine the mapping between Hebbian learning coefficients (weights) to coupling resistors, knowing that they are essential elements for the network functionality. In this work, we proposed a mapping function that translates Hebbian signed weights to coupling resistances in a VO 2 -based ONN for systematic ONN analysis. Our simulations on 60-oscillators test case study highlighted a strong dependency between the ONN recognition accuracy and its coupling resistance range, set by mapping parameter β. Although we identified a suitable value β = N/32 to achieve good ONN accuracy, our mapping formulation provides coupling resistances that differ only with few percent. This would lead to significant FIGURE 14 | (A) Two oscillators coupled by R C , which can be decomposed in two series resistances: R C = R min C + R. (B) Evolution of R, R min C with respect to β. R max ≈ 10% R min C gives the best accuracy results. Note that R min = 1.7% R min C . (C) ONN recognition accuracy for different values of β. The dashed line is obtained for a linear fit of the mapping µ N , i.e., with coupling resistances that are linearly spaced. (D) Impact of R C 's mean variation on recognition accuracy. hardware design constraints, as resistor mismatches smaller than 1.7% would be required to emulate two consecutive weights. In our mapping formalism, we used the phase transition function ζ , which provides the coupling resistance range holding the ONN memory. As we only derived ζ from oscillator dynamics, we believe the oscillator design could be optimized to expand the coupling resistance range and relax synaptic design constraints. For example, oscillator biasing current and supply voltage are the knobs that could be adjusted to maximize the synaptic range.
Here, we reported on a mapping function to compute coupling resistances from signed Hebbian coefficients in a VO 2 -based ONN. We first enhanced the ONN initialization control based on a simple architecture where every oscillator can be decoupled from the network via a switch. We were able to derive the phase transition function from the ONN dynamics, which is crucial for ONN memory. We then merged this analytical formulation with a sigmoid activation function from ANN formalism to build a mapping function. To demonstrate the ONN architecture's applicability with the proposed mapping function, we presented a test case of pattern recognition with 60 fully coupled oscillators. Finally, we showed that ONN recognition accuracy is very sensitive to relative variations between coupling resistances.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.