Distributed coupling complexity in a weakly coupled oscillatory network with associative properties

We present a novel architecture of an oscillatory neural network capable of performing pattern recognition tasks. Two established strategies for obtaining associative properties in oscillatory networks invoke either a physical, time constant or a global, dynamical all-to-all coupling. Our network distributes the complexity of the coupling between the spatial and the temporal domain. Instead of O(N2)?> physical connections or a global connection with O(N2)?> frequency components, each of the N oscillators receives an individual coupling signal which is composed of N − 1 frequency components. We demonstrate that such a network can be built with analog electronic oscillators and possesses reliable pattern recognition properties. Theoretical analysis shows that the scalability is in fact superior to the dynamic global coupling approach, while its physical complexity is greatly reduced compared to the individual time constant coupling.


Introduction
One outstanding property of biological neural networks is the ability to perform pattern recognition tasks. A well known mathematical model emulating the behavior of neurons on an abstract level is the Hopfield network [1] which consists of N neurons each connected to all other neurons individually via interconnections with synaptic strength s i j . Also a network of coupled Kuramoto oscillators [2] may exhibit associative properties similar to the Hopfield network [3][4][5][6][7][8][9][10][11][12]. In particular, it was shown that a system of identical oscillators can recognize binary patterns [8,10]. In a frame of reference rotating with the natural frequency of the oscillators, the evolution equation of each phase variable ϑ i of the oscillatory network readsθ i = N j=1 s i j sin(ϑ j − ϑ i ). (1) A severe obstacle to any hardware implementation of such a network is that at least N (N − 1) experimental parameters representing the coupling matrix s i j are required. Hoppensteadt and Izhikevich [10] suggest an alternative approach based on oscillators with distinct frequencies and a global time dependent coupling of all neural oscillators. Applying a single external coupling function to all oscillators replaces a hard-wired all-to-all coupling by a dynamical all-to-all coupling. Hence, the number of coupling parameters is reduced to just one, which now is time dependent. This concept was validated experimentally with simple electronic oscillating circuits [13]. However, the approach has the disadvantage that the N (N − 1)/2 frequencies present in the global coupling function must be arranged according to a so-called Golomb ruler [14,15], which again restricts the scalability to large systems severely. The needed frequency range scales with N 2 : ω ω < 1 N 2 , with ω being a typical oscillator frequency and ω its accuracy [13].
In this work, we also consider connected neurons as weakly coupled oscillators, but we use a qualitatively different network architecture. Our network invokes N time-dependent coupling functions being composed of N − 1 frequency components each. In comparison, the hard-wired all-to-all coupling uses N 2 coupling functions which are constant in time, while the dynamical all-to-all coupling uses just a single time-dependent coupling function with N (N − 1)/2 frequency components. Below, we will demonstrate that our proposed architecture is more feasible for implementation as hardware than both alternative networks discussed above. In particular, we will quantify that the restrictions on the oscillator frequencies are less severe than for the dynamically coupled network, while the number of coupling parameters scales only linearly with the number of oscillators. Furthermore, we will present a robust experimental implementation of this concept with electronic oscillators. Figure 1 shows the general idea: N oscillators with different natural frequencies ω i and phases ϑ i are subject to an individual coupling signal, which is generated by a central unit. This unit receives the phase information of all N oscillators in the form of their phase shifts ϕ i = ϑ i − ω i t and generates the N coupling functions based on this information. Later, we will show that this system is also effectively governed by equation (1) if the phase ϑ i is replaced by the phase shift ϕ i . System (1) has fixed points (ϕ i − ϕ j ) * = 0 or π depending on the sign of s i j , with s i j = ξ i ξ j and ξ i , ξ j = ±1. Therefore, ξ i and ξ j can be regarded as pixels of a binary pattern ξ with, e.g. ξ i = 1 being black and ξ i = −1 being white. This results in the differences of the phase shifts ϕ j − ϕ i of 0 for the same color (ξ i = ξ j ) and π for the complementary one (ξ i = −ξ j ). For the pattern recognition process we provide a selection of M patterns ξ k by applying the Hebbian rule: The process of pattern recognition is the same as the one used in [10,13].
Note that the coupling strategy depicted in figure 1 has been mentioned before in conjunction with linearly increasing oscillator frequencies ω 1 , ω 1 + , ω 1 + 2 , . . . , ω 1 + N [10]. It turns out, however, that the coupling with this distribution of frequencies is not feasible, which will be elaborated below.

Theoretical considerations
We will now demonstrate that the dynamics of our coupling strategy is also governed by equation (1). Having a network of electronic oscillators in mind, we start with the following system of differential equations for each oscillator: Here, U i is the voltage of oscillator i, I i is the current flowing through the inductance of the van der Pol oscillator (see appendix A), and a i (t) is the individual coupling function. f (U i , I i ) and g(U i , I i ) describe the uncoupled oscillator and ε is a small parameter. Since the coupling is weak, we can describe the evolution of each oscillator with a single differential equation using the concept of the infinitesimal Phase Response Curve (iPRC) [16]: We now proceed by averaging equation (4) over time as we are only interested in the long-term behavior that varies significantly: We use and assume the shapes of the voltages of each oscillator U i (t) to be sinusoidal and ∂ϑ i ∂U i to be cosine shaped. See appendix B for a measurement of the iPRC of the oscillators which justifies this assumption. Note that in equation (6), for each a i (t) the summation runs only over one index ( j), as opposed to the formulation of the global coupling function a(t) [10,13], where the summation runs over both, i and j. Assuming that all oscillators have identical amplitudes, we can derive restrictions on the frequencies of the oscillators that, if fulfilled, transform equations (5) intȯ with c being a constant of magnitude 1. See appendix C for the detailed calculation. Equation (7) matches the desired equation (1). The restrictions on the oscillator frequencies are The first requirement simply states that all frequencies have to be different. The second one is fulfilled if the smallest frequency in the system is greater than one third of the largest frequency. The third one determines the scalability of the system. The set of integers that fulfill the third requirement is composed of all integers that do not contain any 2 in base-3 notation [17], which is de facto a Cantor set. One can show that the space required by the frequencies grows with a power of ln 3 ln 2 : Comparing these requirements with the ones obtained for the single coupling function network [10], one can see that they differ only in the third one. For a purely global coupling all frequency differences need to be different, which is the definition of a Golomb ruler: As already mentioned, this leads to a scaling of the frequency range with a power of 2. Thus, with N individual coupling functions instead of one global coupling function, we gain a better scalability of the system. Moreover, each individual coupling function contains significantly less frequency components and is therefore easier to create.
To get a better intuition for the improved scalability, suppose that you had oscillators with an accuracy sufficient for a network of N global = 10 globally coupled oscillators. Then, you would be able to connect N dist = 18 oscillators with the distributed coupling at hand. For N global = 100, we have N dist = 334, and finally for N global = 1000, the network with the distributed coupling could be as large as N dist = 6103. Hence, with the same frequency accuracy, much larger patterns can be processed in the network with distributed coupling.

Experimental setup
A schematic of the electric circuitry used for the experiments can be seen in figure 2. We used eight van der Pol oscillators, which consisted of a capacitance, a gyrator circuitry instead of a conventional inductor and a negative impedance combined with two diodes instead of a tunnel diode. See appendix A for the complete circuitry and the resulting differential equations. The summation over the voltages of the oscillators was done via an adder circuit realized with an operational amplifier. The added signal was multiplied with the individual coupling functions a i (t) = j =i s i j cos(ω j t − ω i t) using an analog multiplier.
In this experimental demonstration, the coupling signals a i (t) were calculated with a PC according to the measured prevailing frequencies of all oscillators and converted to analog signals to be passed to the multiplier. The multiplied signals were then fed back to the individual oscillators as a small current, dimensioned by the 100 k resistance. Thus, in our realization, the coupling functions a i (t) were generated digitally by a PC, whereas all other computational steps were implemented with analog circuitry.
The experiments were conducted as follows. At first, we determined the effective frequencies of the oscillators by counting the zero crossings of the voltage signal during a defined time interval and dividing the number by the length of this interval. Based on these results, the coupling functions a i (t) were computed. Immediately afterwards the pattern recognition experiment was carried out. It consisted of two parts, the initialization and the recognition part. During the first half of the measuring interval, the pattern ξ which was to be recognized, was initialized. Therefore, the coupling functions a i (t) = 8 j =i ξ i ξ j cos(ω j t − ω i t) were used. As a consequence, all phase shifts ϕ i adopted to one of two values that were π apart, representing the pattern under investigation. One of these values was arbitrarily set to zero. During the second fourth of the measuring interval, we again determined the frequencies as during this time interval the phase shifts do not change anymore. These frequencies were used to obtain the evolution of phase shifts over the entire measuring interval by comparing the actual zero crossings of the voltage signals with the expected zero crossings.
During the second half of the measurement, the recognition was performed. The coupling was switched to the Hebbian superposition of the M patterns ξ k , given by s i j = M k=1 ξ k i ξ k j . The expected result was that the system evolves from the initialized state ξ to a state corresponding to the pattern ξ k most similar to ξ . During this process, the oscillators corresponding to the defective pixels changed their phase shift by π .
All patterns used for recognition were orthogonal, i.e.   Figure 3 shows two exemplary measurements, a typical one (figure 3(a)) and the least favorable one ( figure 3(b)) of 100 experiments performed with randomized patterns and frequencies between 25 and 70 kHz that were chosen from a Cantor set to fulfill the restrictions derived above, namely ω 1 = 25 kHz, ω 2 = 25 kHz + , ω 3 = 25 kHz + 3 , ω 4 = 25 kHz + 4 , ω 5 = 25 kHz + 9 , ω 6 = 25 kHz + 10 , ω 7 = 25 kHz + 12 , ω 8 = 25 kHz + 13 = 70 kHz with = 3461 Hz. During the course of the experiments, frequencies would deviate from these values up to ±100 Hz. The initial pattern ξ (shown in the top left corner) was compared to the three patterns ξ 1 , ξ 2 and ξ 3 (shown in the top right corner). Least favorable here means that the average deviation of all oscillators from the hypothetical ideal behavior (i.e. the defective oscillator changes its phase shift by π , all others do not change their phase shift) was the largest for this experimental run.

Experiments with oscillator frequencies chosen from a Cantor set
In figure 3(a), pattern ξ differs from pattern ξ 1 in the sixth pixel, which is attributed to the magenta colored oscillator. After the change of the coupling at t = 0, this oscillator switches from the '0-branch' to the 'π -branch', whereas the other oscillators stay on their branches: pattern ξ 1 is recognized. We can also observe a constant upward drift, which we do not yet understand completely. However, as one can see, this does not affect the process of pattern The oscillators that are supposed to change are colored in red (hatched), the others in black (not filled). The height of the black columns is given by the axis on the left, the height of the red columns is given by the axis on the right, which is scaled down by a factor 7 (as one out of eight oscillators is defective for each run). The outliers originate from the measurement shown in figure 3(b).
recognition as it operates on a much slower time scale than the phase shift change and, moreover, it is constant and of equal strength for all oscillators. Figure 3(b) shows the worst out of 100 experiments. Yet, one can see that the pattern is still recognized, as there are two clearly distinguishable branches and the correct (orange colored) oscillator switches its branch. Figure 4 shows the histogram of the changes in phase shift ϕ for the series of 100 measurements with randomized orthogonal patterns, which also contained the measurements of figure 3. The phase shifts were evaluated 0.04 s after the coupling was switched to recognition mode.
The resulting distributions for the oscillators that should switch branches and for the oscillators that should remain on their branches center close to the expected values ϕ = π and 0, respectively, and are well separated. The small shift of the left peak toward higher values of ϕ results from the already mentioned upward drift.
For the other peak, this effect is less pronounced because ϕ is slightly smaller than π. This happens because the seven remaining oscillators experience a small recoil effect due to the branch switching of the oscillator corresponding to the defective pixel.
Independent of these details, the important result of this series of measurements is that every single pattern recognition measurement out of 100 was successful.

Experiments with evenly distributed oscillator frequencies
Having shown that a network with frequencies distributed according to equation (8) may recognize patterns reliably, it would be worthwhile to also validate the scaling law equation (9) experimentally. Given the considerable effort this would require, it is more practical to demonstrate that the network fails when the frequencies are chosen in a way that violates equation (8). This is especially interesting since inevitable inaccuracies in the experiment might reduce the effect of undesired resonances. In the following, we show that it is not sufficient to distribute the frequencies evenly ω 1 , ω 1 + , ω 1 + 2 , . . . , ω 1 + N as an intuitive choice that has also been mentioned before [10]. Figure 5 shows two exemplary pattern recognition experiments with evenly distributed frequencies between 25 and 70 kHz, with deviations of up to ±50 Hz. Except for the evenly spaced frequencies, the experiments were conducted in the same way as those depicted in figure 3. Obviously, the oscillators do not settle for a constant distribution of phase shifts during initialization, let alone recognition 2 . Instead, the time evolution of the phase shifts fluctuates with amplitudes of up to about π/2 due to spurious resonances. Figure 5(a) shows a successful pattern recognition as the oscillator corresponding to the red curve is the only one exhibiting a substantial change in phase shift ϕ after the change of coupling. In figure 5(b), on the other hand, the oscillator corresponding to the defective pixel (orange) does not exhibit a substantial change in phase shift, while others do (red, green, cyan). The oscillators that are supposed to change are colored in red (hatched), the others in black (not filled). The height of the black columns is given by the axis on the left, the height of the red columns is given by the axis on the right, which is scaled down by a factor 7 (as one out of eight oscillators is defective for each run). Figure 6 is comparable to figure 4. It shows the histogram of phase shift changes for a series of 100 measurements with randomized patterns and evenly distributed frequencies, which include the measurements of figures 5(a) and (b). In contrast to figure 4, the two distributions overlap considerably, which results in the failure of more than 40 out of the 100 pattern recognition runs. This serves to show that the frequency conditions derived in section 2 are indeed relevant for the experimental system at hand, even in the presence of inherent frequency inaccuracies, which tend to alleviate spurious resonances.

Conclusion
In this paper, we have presented a new type of weakly coupled oscillatory network which is capable of recognizing binary patterns. The network allows for larger system sizes than the network discussed in the literature [10,13]. In addition, it is more feasible to implement as hardware than traditional artificial neural network architectures. A simple experimental setup was introduced which is capable of performing pattern recognition tasks reliably. The results are a promising foundation for further development, involving, e.g. more sophisticated oscillators or the on-chip generation of the coupling function. They also point to interesting fundamental questions concerning the dynamics of oscillatory networks.

Appendix A. Circuit details
In this work we used electronic oscillators which are essentially van der Pol oscillators. We replaced some classic elements of the circuit. In particular, we exchanged the inductor for a so-called gyrator circuit with a capacitor, and the tunnel diode was substituted by a negative impedance with two diodes. Figure A.1 shows the electric circuit of a single oscillator. In the middle of the figure (colored in black) one can see the 1 nF capacitor of the oscillator. To the left (colored in red) is the negative impedance circuitry with the two diodes. This circuit effectively transforms the 6.19 k resistor into a −6.19 k one. Combined with the two diodes, this part in total provides the typical U-I-characteristic of a tunnel diode. To the right (colored in blue) one can see the gyrator which transforms the 33 nF capacitance into an inductance. The frequency adjustments are done via the 1 k potentiometer R var . Consequently, the inductance ranges between L i = 3.3-36.3 mH; this results in a frequency range of 25 to 70 kHz. U i is the voltage fed into the adder circuit (see figure 2, which shows how the coupling of the oscillators is implemented).
The whole circuit of oscillators and coupling is described by the following differential equations:U with C = 1 nF and I negImp being the current flowing through the negative impedance circuit. R and R are the resistances shown in figure 1, with R = 8.25 k and R = 100 k . The denominator in the coupling term contains a factor 10, as the AD633 Multiplier automatically divides by 10.
If we compare equations (2) and (3)   For the measurement of the PRC we provided a short high amplitude voltage pulse (see equation (B.2) below) to one otherwise undisturbed van der Pol oscillator and determined the change in phase 294 times at random oscillator phases.
We connected the node at potential U i in figure A.1 to an external voltage source U ext which provides the pulses via a resistance R = 8.25 k . Therefore, the voltage equation for the oscillator readsU (with C = 1 nF, R = 8.25 k and R = 100 k ). The amplitude of the oscillator was 0.52 V with a frequency of 61.5 kHz. The voltage pulse we used had the shape (B. 2) The expected voltage jump induced by this pulse is U = U ext RC dt ≈ 0.02 V. To determine the phase response Z ϑ = ϑ U we calculated the change in phase ϑ by comparing the time between two zero crossings of the voltage before and right after the pulse.
Further experimental details are identical to the ones in reference [13]. Figure A.1 shows that it is justified to model the PRC by a cosine function.

Appendix C. Averaging the phase equations
In the following, we show how averaging reduces equations (4)- (7). We start with the following equation:θ Here, ∂ϑ i ∂U i = Z max PRC,i cos ϑ i is the PRC, with Z max PRC,i being the amplitude and U j = U max j sinϑ j is the voltage of oscillator j with U max j being its amplitude. Using ϑ = ωt + ϕ, equation (C.1) transforms tȯ · · · + sin((ω j + ω i )t + ϕ i + ϕ j )) dt. (C.2) Now we insert the coupling function a i (t) and transform the equation via trigonometric identities: Since all frequencies are positive, the integral over the third term is always 0. The integral over the first term is 0 if ω i = 1 2 (ω k + ω j ) is fulfilled. Furthermore, the fourth term vanishes if ω i = 1 2 (ω k − ω j ). This leaves only the second term. If we guarantee all frequencies to be different, we obtain a non-vanishing contribution if k = j. This results iṅ If we assume all the voltage amplitudes U max j to be equal, we obtain the desired form. The constant c in equation (7) is thus c = 1 4 Z max PRC,i U max . Ideally, Z max PRC,i U max = 1 should hold. During our measurement of the PRC (see appendix B) we found U max = 0.52 V and Z max PRC = 1.77 rad V −1 (curve fit).