A Complex-Valued Oscillatory Neural Network for Storage and Retrieval of Multidimensional Aperiodic Signals

Recurrent neural networks with associative memory properties are typically based on fixed-point dynamics, which is fundamentally distinct from the oscillatory dynamics of the brain. There have been proposals for oscillatory associative memories, but here too, in the majority of cases, only binary patterns are stored as oscillatory states in the network. Oscillatory neural network models typically operate at a single/common frequency. At multiple frequencies, even a pair of oscillators with real coupling exhibits rich dynamics of Arnold tongues, not easily harnessed to achieve reliable memory storage and retrieval. Since real brain dynamics comprises of a wide range of spectral components, there is a need for oscillatory neural network models that operate at multiple frequencies. We propose an oscillatory neural network that can model multiple time series simultaneously by performing a Fourier-like decomposition of the signals. We show that these enhanced properties of a network of Hopf oscillators become possible by operating in the complex-variable domain. In this model, the single neural oscillator is modeled as a Hopf oscillator, with adaptive frequency and dynamics described over the complex domain. We propose a novel form of coupling, dubbed “power coupling,” between complex Hopf oscillators. With power coupling, expressed naturally only in the complex-variable domain, it is possible to achieve stable (normalized) phase relationships in a network of multifrequency oscillators. Network connections are trained either by Hebb-like learning or by delta rule, adapted to the complex domain. The network is capable of modeling N-channel electroencephalogram time series with high accuracy and shows the potential as an effective model of large-scale brain dynamics.


INTRODUCTION
Currently, there are two prominent approaches to characterizing the neural code: the instantaneous rate or frequency at which the neuron fires action potential ("the spike frequency code" or the "rate code") and the time of the occurrence of action potential ("spike time code"). The former assumes the information is processed/encoded over a larger time scale conveyed by average number of spikes fired in a given duration forms the basis of a large class of neural networks called the rate coded neural networks (Lippmann, 1990;Ruck et al., 1990;Lawrence et al., 1997). Whereas the latter believe that the precise timing of the action potential fired by a neuron encodes the ongoing activity, giving rise to a broad class of spiking neural network (Maass, 1997a,b;Izhikevich, 2003Izhikevich, , 2004Ghosh-dastidar and Lichtenstein, 2009). These two classes of neural networks are capable of universal function approximation and well explored (Maass, 1997a;Auer et al., 2008). They have been reported to solve a wide range of information processing problems like vector space transformation, dimensionality reduction, sequence processing, memory storage as an attractor and autoencoding (Hopfield, 1982;Frasconi et al., 1995;Kohonen, 1998;Trappenberg, 2003;Schmidhuber, 2015).
The third class of neural code, dubbed the oscillation coding (Fujii et al., 1996) also has been around for more than two decades. However, it has not been adopted in modeling literature as extensively as the first two types of neural code. Oscillation coding is based on the idea that not single-neuron activity, but the synchronized collective activity of a neural ensemble is the true building block of the brain's activity. Unlike the rate code, which is often represented as a continuous signal encoding information as an average neural activity, or the spike time code, which is often described as a train of delta functions, oscillation code is a smooth signal that is said to be composed of distinct frequency bands. Oscillation coding also offers an opportunity to represent important temporal phenomena like temporal binding via transient synchrony Terman, 1995, 1997;Shi et al., 2008) or rhythmic behaviors such as locomotion (Ermentrout and Kopell, 1989). There have been comprehensive and valorous attempts to describe all brain function in terms of the oscillatory activity of neural ensembles (Draguhn, 2004;Buzsáki et al., 2012). However, since existing literature does not commit to the exact size of a "neural ensemble, " activity measured at different scales goes by different names, including local field potentials (LFPs), electrocorticograms (ECoGs), electroencephalograms (EEGs), etc.
Since oscillations are ubiquitous in the brain, one would naturally expect nonlinear oscillator models to be used extensively to describe brain function. However, models of brain function often use sigmoidal neurons (to represent the rate code), spiking neuron models (to represent the spike time code), or more detailed conductance-based or biophysical models. In the field of neural signal processing, for example, a large body of literature depicting various deep neural network architectures such as recurrent neural network, convolutional neural network and deep belief network have been used for EEG signal processing and classification at various scales which found application in the following areas: brain-computer interface, epilepsy, cognitive and effective monitoring, etc. (Craik and He, 2019;Roy et al., 2019). Nonlinear oscillators are used exclusively for describing oscillatory phenomena like rhythmic movements or feature binding by synchronization. In an attempt to correct this bias, it would be a worthwhile modeling exercise to construct general networks of nonlinear neural oscillators that can be used to describe a wide range of brain functions.
The coupling strategies between a pair of neurons differ depending on the neuronal dynamics defined by the mathematical representation of the neural code. Biophysical neuron models employ biophysically detailed synapse models for coupling whereas the spiking neuron models employ elementary single exponential alpha function or double exponential function for coupling (Destexhe et al., 1994). Various two-variable oscillator models such as FitzHugh Nagumo model, Van der Pol, Wilson-Cowan, Kuramoto and Hopf oscillator principally employ what may be described as "real coupling strategy" (Campbell and Wang, 1994;Toral et al., 2003;Low et al., 2005;Izhikevich, 2007;Breakspear et al., 2010;Hoff et al., 2014;Sadilek and Thurner, 2015). In real coupling strategy, only the main variables (x 1 and x 2 ) of the two oscillators interact directly. Delayed real coupling strategies have also been explored for several of these oscillator models (Wirkus and Rand, 2002;Saha and Feudel, 2017). The averaged behavior of pulse-coupled oscillators has also been analyzed (Ermentrout and Kopell, 1991). On the contrary, in the complex coupling strategy (Hoppensteadt and Izhikevich, 2000), the two variables of the oscillator (x and y) are combined to form a complex variable (z = x + iy), and the complex state of one oscillator influences the complex state of the other. In this case, both x and y variables of one oscillator influence the x and y variables of the other.
In this paper, we describe a general trainable network of neural oscillators. For the oscillator model, we choose the simple Hopf oscillator. One reason behind the choice of the Hopf oscillator is the elegant form its description assumes in the complex variable domain. We show that it is important to operate in the complex domain in order to define some of the novel modeling features we introduce in this paper. Particularly, we introduce the concept of "power coupling" by which it is possible to achieve robust phase relationships among oscillators with very different intrinsic frequencies. The proposed network of oscillators can be trained to model multi-channel EEG, thereby demonstrating the potential to evolve into a large-scale brain model.
The outline of the paper is as follows. Section "Prior Work: Complex Coupling" begins with the definition of Hopf oscillator in the complex domain and show how to adapt the intrinsic frequency of the oscillator to that of a sinusoidal forcing signal. Section "A Pair of Hopf Oscillators With Complex Coupling" describes the dynamics of a pair of coupled Hopf oscillators and shows the advantages in coupling with a complex coupling constant. Section "Coupling Two Hopf Oscillators With Different Natural Frequencies" highlights the difficulties in coupling a pair of Hopf oscillators with very different intrinsic frequencies, and Section "Power Coupling" Between a Pair of Hopf Oscillators" shows how the difficulties can be overcome by adopting "power coupling." Gathering all the modeling elements developed so far, Section "Adaptive Hopf Oscillator With Complex Input or Complex Adaptive Hopf Oscillator" describes a network of oscillators that can model multi-channel EEG. A discussion of the work was presented in the last section.
representation state variable : z = re i∅ ; ; i = √ −1 , Cartesian coordinate representation state variables : x, y and polar coordinate representation (state variables : r, ∅) respectively: Complex state variable representation: Cartesian coordinate representation: Polar coordinate representation: Depending on the parameter values (α, β 1 , β 2 ), the autonomous behavior of the oscillator principally falls into four regimes: critical Hopf regime (α = 0, β 1 < 0, β 2 = 0), supercritical Hopf regime (α > 0, β 1 < 0, β 2 = 0), supercritical double limit cycle regime α < 0, β 1 > 0, β 2 < 0, local maxima > 0 and subcritical double limit cycle regime α < 0, β 1 > 0, β 2 < 0, local maxima < 0 (Kim and Large, 2015). The critical Hopf regime has a stable fixed point at the origin and has the ability to show a stable resonating response to the complex sinusoidal input (Fe iω 0 t ). The supercritical Hopf regime has an unstable fixed point at the origin and a stable limit cycle at r = α |β 1 | . In the supercritical double limit cycle regime, the system exhibits two limit cycles, one of which is stable while the other being unstable. In the subcritical Hopf regime, the system has one stable fixed point at the origin. However, it has the ability to show stable oscillation under the influence of complex sinusoidal input whose frequency is not too different from that of the oscillator (ω − ω 0 < ). Throughout the rest of the paper, we will be using supercritical Hopf regime (with α = µ, β 1 = 1, β 2 = 0), which can be defined as follows: Complex state variable representation: Cartesian coordinate representation: Polar coordinate representation: We now present a series of results related to a single oscillator, a coupled pair, and a network of Hopf oscillators in the supercritical regime defined above.

Prior Work: Complex Coupling
Single Oscillator With Adaptive Frequency It has been previously shown that when a Hopf oscillator is influenced by a real sinusoidal input signal, it can adapt its natural frequency to the frequency of the input signal if it follows the following dynamics (Righetti et al., 2005).
Given, z = re ∅ = x + iy. It can be recognized that the Hopf oscillator is perturbed by a real input signal I ext (t) with coupling strength ε. At steady-state r reaches √ µ for ε = 0. If I ext (t) = I 0 sin (ω 0 t + ϕ) (I 0 , ω 0 , and ϕ are the magnitude, frequency and phase offset of the sinusoidal input signal), the natural frequency, ω, adapts to the frequency ω 0 .

Adaptive Hopf Oscillator With Complex Input or Complex Adaptive Hopf Oscillator
When a Hopf oscillator is influenced by a complex sinusoidal input signal, the natural frequency of the oscillator can adapt to the frequency of input if the natural frequency of the oscillator is updated according to the following equation: In this scenario, I ext (t) is a complex sinusoidal signal. It is straight forward to derive the learning rule for the natural frequency of the oscillator if eq. 2a is represented in the Cartesian and polar coordinate forms, respectively, as follows:ẋ In the phase plane representation, it can be observed from eq. 2a2 that the influence caused by the input perturbation on the oscillator phase is εI 0 r sin (ω 0 t − ∅). Whereas from the Cartesian coordinate system representation (eq. 2a1) it can be observed that the overall influence of the external perturbation on the phase point along the tangential axis of the limit cycle ( − → e ∅ ), − → P ∅ is (to understand the vector notations refer to Figure 1): This motivates us to adapt the learning rule for the natural frequency of the oscillator, as proposed in eq. 2c, by dropping the magnitude of oscillation because of the same reason as described by Righetti et al. (2006). We have simulated the eqs. 2a-2c and observed that the proposed learning rule for the natural frequency of the oscillator allows it to adapt to the frequency of the input complex sinusoidal signal, as shown in the following Figure 2.

A Pair of Hopf Oscillators Coupled Through Real Coupling
When two Hopf oscillators with equal natural frequencies are coupled with the real coupling coefficient as described below (eq. 3), they are going to phase lock in phase (0) or out of phase by (2n + 1)π depending on the polarity of the coupling and initialization. (for proof refer to Supplementary Appendix 1) where W 12 is the real coupling coefficient from 2nd oscillator to 1st one, and W 21 is the real coupling coefficient from 1st to 2nd. In the stated scenario, both of the oscillators have identical natural frequencies, ω 1 = ω 2 . This, with real-valued coupling, a pair of Hopf oscillators with equal intrinsic frequencies can only produce two possible values of phase difference.

A Pair of Hopf Oscillators With Complex Coupling
When two Hopf oscillators with identical natural frequencies are coupled bilaterally through complex coefficients with Hermitian symmetry, they can exhibit phase-locked oscillation at a particular angle similar to the angle of complex coupling FIGURE 1 | Considering the circle as the limit cycle of the Hopf oscillator, e ∅ is the unit vector along increasing azimuth angle (∅), e r is the unit vector along radius (r), P is the overall external input perturbation, P ∅ and P r are the external input perturbation along e ∅ and e r respectively. The analogy is drawn from Righetti et al. (2006).

Simulated Result
To verify the above result, we have simulated eqs. 4a1 and 4b1 numerically by Euler integration method on MATLAB platform with t = 1 msec. It can be observed from the plot (refer to the Figure 3) that the steady-state phase difference between the two oscillators (ψ ss ) is either achieving θ or 2π + θ depending on the whether it was initialized in the intervals θ − π < ψ (0) ≤ θ + π and θ + π < ψ (0) ≤ 2π + θ respectively.

Training Rule for the Complex Coupling
When a pair of complex sinusoidal inputs with identical frequency is presented to a pair of coupled Hopf oscillators with complex coupling and identical natural frequencies (eq. 4), and when the complex coupling coefficient is adapted according to eq. 5, the angle of the coupling coefficient approaches the phase difference of the external inputs. In other words, the complex coupling learns the phase relationship of the external inputs.
On the other hand, in a network of two coupled oscillators driven by separate external sinusoidal forcing as described in the following eq. 6 with the same frequency as the natural frequency of the two oscillators but any phase offset (ϕ 1 , ϕ 2 ) will drive those oscillators to oscillate at the same phase as the external sinusoidal forcing (∅ 1 ≈ ωt+ϕ 1 and ∅ 2 ≈ ωt+ϕ 2 ), provided very low magnitude of the complex coupling coefficient (A 1).
Under the influence of external input, the angles of z 1 and z 2 (∅ 1 and ∅ 2 respectively) tend to ωt + ϕ 1 and ωt + ϕ 2 , respectively. When the complex coupling coefficient is trained according to eq. 5, the angle(W) or θ learns angle (z 1 z 2 * ) or ∅ 1 (t) − ∅ 2 (t) ≈ ϕ 1 − ϕ 2 . During training, A must be kept low (A 1) so that z 1 and z 2 dynamics is influenced only by the intrinsic oscillatory dynamics and the forcing signal. During testing, however, the external input is removed, A is again increased so that the phase difference between the oscillators is stable and determined by the coupling constant.

Simulated Result
The proposed Hebb like learning rule (eq. 5) allows the angle of complex coupling weight (θ) to learn the phase difference between the complex sinusoidal inputs driving each of the oscillators for very low magnitude of the complex coupling weight (A 1). To verify this, eq. 6 is simulated for the set of parameters as described in the following Figure 4. It can be observed that θ learns the phase difference between the complex sinusoidal input signals (ϕ 1 − ϕ 2 ).

Coupling Two Hopf Oscillators With Different Natural Frequencies
A different kind of dynamics can be observed when two sinusoidal oscillators with different natural frequencies are coupled with real coupling coefficient. A pair of Hopf oscillators with real coupling shows Arnold tongue behavior which may be described as follows: the two oscillators entrain to commensurate frequencies (having a simple integral ratio) for specific ranges of the coupling strength, and the range of frequencies for which the entrainment occurs widens with the strength of the coupling coefficient (Izhikevich, 2007). There are whole ranges of the coupling strength and the natural frequencies where there is no specific phase relationship at all. Therefore, it is not easy to get a stable phase relationship between two Hopf oscillators with very different frequencies and real coupling.
The situation is a bit more facile with Kuramoto oscillators (Izhikevich, 2007). When a pair of Kuramoto oscillators are coupled according to the equation shown below, their natural frequencies converge to a particular frequency in between their FIGURE 4 | (A) Network schematic of dynamics represented in eq. 6 where the complex coupling coefficient is trained according to eq. 5. (B) It can be verified that for the following set of network parameters θ learns ϕ 1 − ϕ 2 .
natural frequencies, depending on the coupling strength.
Let, K 1 and K 2 be the real coupling strengths in the forward and backward directions, respectively. The phase difference between the two oscillators (∅ = ∅ 1 − ∅ 2 ) reaches a steady-state value if K 1 + K 2 > |ω 1 − ω 2 | which is satisfied by the relationship sin ∅ * = ω 1 −ω 2 K 1 +K 2 where ∅ * be the steadystate phase difference. Notably, both of the oscillators reach a periodic solution with frequency ω * = K 1 ω 2 +K 2 ω 1 K 1 +K 2 . In this scenario, although synchronization at a particular intermediate frequency with a particular phase difference is ensured, the natural frequency of oscillation is not maintained, i.e., to phase lock at a particular angle, the oscillation of both the oscillators had to converge at a particular intermediate value. A network of such Kuramoto oscillators (eq. 7a) with natural frequencies drawn from a Gaussian distribution with constrained standard deviation can show increased phase synchrony among the oscillators when the isotopic coupling strength (K) exceeds a threshold (Breakspear et al., 2010).
Thus, we can see that the dynamics of a pair of coupled Hopf oscillators become dramatically more complicated when we relax the equality relationship between the natural frequencies of the coupled oscillators. We can visualize a simple, desirable extension of the dynamics of a pair of coupled oscillators with distinct natural frequencies. When the natural frequencies are equal, we found that the phase difference equals the angle of the complex coupling factor. However, when the natural frequencies are unequal, there cannot be a phase difference in the usual sense (∅ 1 − ∅ 2 ), although when the coupled oscillators are entrained in the m : n ratio, the phase difference has been defined as m∅ 1 − n∅ 2 (where, m, n are integers). Therefore, one may define a quantity called "normalized phase difference" as follows: Let ψ 12 and ψ 21 be the normalized phase difference of 1st oscillator w.r.t 2nd oscillator and vice versa. Is it possible to relate the normalized phase difference with the angle of the complex coupling coefficient? It turns out that such a relationship is not possible even with the complex coupling of eq. 4. To this end, we extend the complex coupling to a new form of coupling we label "power coupling" as described below.

"Power Coupling" Between a Pair of Hopf Oscillators
A pair of sinusoidal oscillators (Hopf or Kuramoto oscillator) can entrain at a specific normalized phase difference if they are coupled through complex "power coupling" coefficient according to the following eqs. 8a and 8b at a particular value dependent on the angle of complex "power coupling" coefficient, natural frequencies of the coupled oscillators and the initial values of the phase angle of those oscillators (for proof please refer to Supplementary Appendix 4).
Coupling a pair of Hopf oscillators through power coupling: ω 2 =W 12 is the weight of power coupling from 2nd oscillator to the 1st oscillator and A 21 e i θ 21 ω 1 =W 21 is the weight of power coupling from 1st oscillator to the 2nd oscillator.
The polar coordinate representation: The schematic of the coupling architecture is elaborated in the following Figure 5. Based on eqs. 8a2 and 8b2, we understand that two modified Kuramoto oscillators may be coupled using power coupling as follows:∅

Power Coupling Among N Hopf Oscillators
Similarly, the dynamics of the i th oscillator in a network of N supercritical Hopf oscillators coupled through power coupling can be represented as: ω j be the weight of power coupling from the j th oscillator to i th oscillator. The polar coordinate representation is: FIGURE 5 | The schematic network representation of eq. 8. There is a kind of frequency transformation occurring at the power coupling synapse because of the following term in the dynamical eq. 8: . At the power coupling synapse, this frequency-transformed version of the oscillation from j th oscillator is weighted by the power coupling weight W ij = Ae i θ ij ω j to convert the forcing signal from j th oscillator to i th oscillator and vice versa.
Frontiers in Computational Neuroscience | www.frontiersin.org Therefore, the dynamics of the i th oscillator in a network of N number Kuramoto oscillators coupled through power coupling can be represented as: We may appreciate the difference in the second term in eqs. 7a and 12.ψ ij s (defined in eq. 7b) can be obtained from eqs. 8a2, 8b2 and 11b for a network of two oscillators and also for the general case of N(> 2) Hopf oscillators respectively: At steady state, ψ 12 and ψ ij = ∅ i ω i − ∅ j ω j can attain any of the possible solutions of the above eqs. 13a,b.
For a network of any arbitrary N number of Hopf oscillators as described in eqs. 10 and 11 under the constraint e iθ ij =e −iθ ji or θ ij =−θ ji ,A ij =A,µ= 1, ψ ij s can achieve the following desirable solution: for certain initial conditions. However, there can exist other solutions for which the sum of the terms on the right-hand side of eqs. 13a,b is zero, without all the individual terms being zero. That is eq. 13b.s above is not satisfied. We call these solutions spurious solutions.
Whether the final solution is a desirable or a spurious solution depends on the initial conditions. Therefore, the solution ψ ij * can be achieved only for certain initial values of ∅ i s.

Three Hopf Oscillators Coupled Through Power Coupling
We now explore the solution space numerically for a threeoscillator system. When three Hopf oscillators are coupled according to the eq. 10 or 11 (N = 3) with θ ij s, chosen as , at steady state σ ij achieves one of the possible solutions depending on the initializations (∅ i (0) s). The∅ i dynamics of the three oscillators can be expressed in the following form.
To verify this numerically eq. 11 was simulated for the given set of network parameters: ω 1 = 4, ω 2 = 7, ω 3 = 10, A ij = A = ∅ i s in the range 0 < ∅ i (0) ≤ 2π with ∅ i (0) = 0.1 for 100 s to let the system achieve a steady-state. It can be observed from eq. 14 that the possible solutions of σ ij ss do not depend on θ ij but on the natural frequencies of the oscillators (ω i ). The following Figure 7 summarizes the dependency of σ ij ss on ∅ i (0) and ∅ j (0) where it can be observed that σ 12ss , σ 23ss , and σ 31ss achieve any of the 5, 3, 4 solutions respectively (mentioned in the figure caption) and particular combinations of these solutions satisfy the above equation. For a given initial value of ∅ 1 , ∅ 2 and ∅ 3 , σ ij ss 's attained one of the solutions and for a confined subspace of ∅ 1 (0), ∅ 2 (0) and ∅ 3 (0) under space 0 ≤ ∅ i (0) < 2π, σ ij ss 's attained the zero or desired solution σ * ij .

Hebbian Learning for the N-Oscillator System With Power Coupling
We may extend the Hebbian learning of complex coupling (eq. 5), to the case of power coupling as follows.
The polar coordinate representation: We have limited the scope of our study only to θ ij dynamics by assumingȦ ij = 0. Let τ W be the time constant of the above learning dynamics. It is shown both analytically and numerically that ω j times angle W ij or θ ij tries to learn ∅ i ω j − ∅ j ω i when the weight of the power coupling is trained according to eq. 15 (for proof, please refer to Supplementary Appendix 5). Similarly, as it was previously shown during the training or phase encoding of the complex coupling coefficient with very small fixed value of A ij (magnitude of complex coupling coefficient), θ ij (angle of complex coupling coefficient) learns ∅ i (0) − ∅ j (0) ≈ ∅ i − ∅ j for arbitrary values of ω (constraining the natural frequencies to be identical). When the power coupling weight is updated according to eq. 15 under the same constraint as A ij 1 andȦ ij = 0, θ ij learns ∅ i (0) ω j − ∅ j (0) ω i when there is no external input.
From eq. 11, it can be interpreted that ∅ i (t) ≈ ω i t + ∅ i (0) for A ij 1. So, when a network of sinusoidal oscillators coupled according to eqs. 10 or 14 with the power coupling weights (W ij s) trained according to eq. 15, In such a network when individual oscillators z i s or ∅ i s are driven by complex sinusoidal forcing (I exti = I 0i e i(ω i t+ϕ i ) for the i th oscillator) and the power coupling weights are trained as stated, θ ij learns ϕ i ω j − ϕ j ω i as with A ij 1,Ȧ ij = 0 and I 0i A ij . ∅ i s are approximately same as ω i t + ϕ i . The dynamics of such a network of supercritical Hopf oscillators is defined below: (16) I exti = I 0i e i(ω i t+ϕ i ) Polar coordinate representation: The dynamics of such a network of Kuramoto oscillators:

Simulated Result
When two Hopf oscillators are coupled through power coupling with very weak coupling coefficient A ij 1 as described above in eq. 16, θ ij learns ∅ i (0) ω j − ∅ j (0) ω i while there is no external input (I exti = 0) and θ ij learns ϕ i ω j − ϕ j ω i when i th oscillator is forced by complex sinusoidal external perturbation. The following simulation results support the stated argument (Figure 8).

A Network of Oscillators With Adaptive Frequency and Trainable Lateral Connections
A pair of coupled adaptive Hopf oscillators driven by distinct complex sinusoidal inputs are capable of adapting their natural frequencies to the frequencies of the complex sinusoidal input signals. The trainable power coupling weight can encode the normalized phase difference between the two complex sinusoidal input signals.
When a pair of Hopf oscillators coupled through trainable power coupling coefficient (eq. 15) are driven by complex sinusoidal inputs (eq. 16), they adapt their natural frequencies according to eqs. 2b or 2c, ω j times the angle of power coupling coefficient W ij (power coupling weight from j th oscillator to i th oscillator) or θ ij can learn the normalized phase difference between the complex sinusoidal input signal. If I exti is the complex sinusoidal input signal driving i th oscillator, then the normalized phase difference among them is: When the following dynamics is simulated for N=2 it can be verified that θ ij learns ϕ j ω 0i − ϕ i ω 0j after ω i , the natural frequency of i th oscillator learns the frequency of the input signal ω 0i (described in Figure 9).

A Network for Reconstructing a Signal by a Fourier-Like Decomposition
In the previous sections, we have described a network of oscillators in which the natural frequencies and lateral connections can be trained. We now add a feature to the network of Section "A Network of Oscillators With Adaptive Frequency and Trainable Lateral Connections" and make it learn an unknown signal by performing a Fourier-like decomposition.
FIGURE 9 | Equations 19a-c are simulated for the following set of parameters:ω 01 = 20, ω 02 = 30, ω 1 (0) = 30, ω 2 (0) = 40, A 12 = A 21 = 0.0001, θ 12 (0) = −1.7884, θ 21 (0) = 1.7884, τ W = 10 3 , ε = 0.9 and I ext1 = e i 20t+ π 4 , I ext2 = e i 30t+ π 6 . It can be verified that θ ij learns ϕ j ω 0i − ϕ i ω 0j after the ω i s learns the corresponding ω 0i s.  Righetti et al. (2005) with additional asymmetric power coupling connections between the oscillators, each of which is driven by the same error signal e(t) between the D(t), teaching time-series signal and P(t), the linear summation of the oscillations of each of the oscillators in the reservoir through the real feed-forward connection weights, elaborated in eq. 20. (B) The network architecture of 2 nd phase of training with the same reservoir of oscillators with tuned natural frequencies and trained power coupling weights in the 1 st phase of training. Y p i is the predicted output signal at the i th output node. Complex feed-forward weights W f ij s connecting j th oscillator to the i th output node, are trained in batch mode according to eq. 23.
adaptable central pattern generator (aCPG) has been proposed by Righetti et al. (2005), the detailed comparison of which is discussed in Section "Discussion." To reconstruct the teaching time-series signal at the output summation node, the error signal that drives each of the oscillators drives the phase offset to the desired phase offset associated with the corresponding frequency present in the teaching time-series signal. The activation and learning dynamics of the network are given below.
Let D (t) be the teaching time-series signal with a finite number of frequencies (D (t) = N i=1 a i cos (ω i t + ϕ i ) where ϕ i be the phase offset associated with i th frequency component), z i = r i e i∅ i , α i s are the real feed-forward weights from i th oscillator to the output summation node, τ W , η ω and η α are the time constant of the learning dynamics for power coupling coefficient, the learning rate of the natural frequency of the oscillators and the real feed-forward weights from oscillators to the output summation node respectively, P (t) is the reconstructed signal at the output summation node. The numerical simulation of the proposed network ( Figure 10A describes the schematic) is elaborated in the following Figure 11. We propose a similar network for Fourier decomposition of complex time series signal with a finite number of frequencies (assuming D (t) = N i=1 a i e i(ω i t+ϕ i ) ), which is an extension of the network proposed in section "A Network of Oscillators With Adaptive Frequency and Trainable Lateral Connections." The proposed network is comprised of a reservoir of complex adaptive Hopf oscillators coupled through trainable power coupling weights, and driven by distinct complex sinusoidal inputs with arbitrarily different frequencies. The network is capable of learning the frequencies and encode the normalized phase relationship among the oscillatory components of the input signal. It is driven by the error signal between the complex teaching signal and the linear summation of complex activations of the Hopf oscillators. The dynamics of the network is described in the following equations.
The schematic of the network is identical to the network described in Figure 10A, and the simulated result for a D (t), which contains three frequency components is as follows (Figure 12).

A Generative Network Which Is Capable of Modeling EEG Signals
The proposed network is capable of modeling an arbitrary number of signals with an overlapping frequency spectrum. It is trained in two phases. In the first phase, a network exactly similar to the one in the previous section (Section "A Network for Reconstructing a Signal by a Fourier-Like Decomposition") is used to encode the constituting frequency components of one of the input signals. During this phase, the natural frequencies (ω i s), the real feed-forward weights (α i s) as well as the power coupling weights W ij s of a network of N Hopf oscillators are trained using the same learning rule described in the previous section (eq. 20). One key difference compared to the previous section is that the teaching signal in the present scenario is aperiodic (frequency spectrum is continuous) and has a finite duration. To overcome this issue, the limited duration teaching EEG signal is presented repeatedly to the network over multiple epochs. This helps the FIGURE 12 | Equations 21a-e are simulated for a teaching signal constituting three frequency components with the parameters D (t) = 2e i 4t+ π 2 + 1.5e i(8t+ π 5 ) + 1.8e i 12t+ π 12 , A ij = 10 −5 , η ω = 0.1, η α = 10 −4 , τ W = 10 4 , ε = 0.5, N = 3, dt = 10 −3 sec. (A) The real part of the network learned signal at the output summation node P (t) and the D (t) signal at various 10 s of time interval. (B) The natural frequencies of the three oscillators learn the three frequency components present in the teaching signal. (C) ω j times the angles of the complex power coupling weights θ ij learn ϕ i ω j − ϕ j ω i . (D) The real feed-forward weights α i 's learn a i 's. network to learn some sort of Fourier decomposition of the teaching signal.
In the second phase of training, the trained reservoir of oscillators tries to reconstruct M number of signals at the corresponding M output nodes assuming these M signals are the outcome of the same underlying process (producing signals with frequencies confined to a certain frequency band) by training the complex feed-forward weights connecting N oscillators of the reservoir to the M output nodes. The schematic of the network architecture ( Figure 10B) and the corresponding learning rules are described below.

First Phase of Learning
During the first phase of learning, an identical network with identical learning dynamics as described by eq. 20 is used. The only fundamental difference between the previous and present scenario is that D (t) is a limited duration aperiodic or quasiperiodic signal compared to the previous case where D (t) was an infinite duration periodic or quasiperiodic signal. i.e., there can be an infinite number of frequency components present in D (t) as the frequency spectrum of an aperiodic signal is continuous (refer to Figure 14A). So, the Fourier decomposition using the proposed network can be accomplished by discretely sampling the continuous frequency spectrum.

Second Phase of Learning
In this phase, the supervised batch mode of learning is used to train the complex feed-forward weights, W f , defined as follows.
Here K ij and ζ ij are respectively the magnitude and angle of the complex feed-forward weights. The derivation for the batch update rule for K ij and ζ ij as described in eq. 22 is given in the Supplementary Appendix 6.
(B) After the second phase of training, the network ( Figure 10B) is able to reconstruct the M = 5 output signals (Y d i for i = 1 to M) with identical frequency components but randomly chosen amplitude (A i ) and phase offset (∅ i ) with high accuracy. Learning parameters: η K = 3 × 10 −5 , η ζ = 10 −6 , no. of epochs for the batch mode of learning = 10000.

Simulation Results
At first, the proposed network is simulated to model an arbitrary number of quasiperiodic signals with identical frequency components. The desired signal (D (t)) consisting of three frequency components (provided ω i ω j is not an integer) with a duration of 20 s is used for learning. The network output P (t) signal (after the first phase of learning) and the respective frequency spectrum is plotted in Figure 13A. After the second phase of training, the desired and the network reconstructed signals are plotted in Figure 13B, and the training parameters are described. One crucial condition has to be met in order to achieve an accurate reconstruction: All the oscillators in the reservoir has to be initialized at z i (0) = 1.
The network performance is tested on low pass filtered (cut off frequency of 5 Hertz) EEG data collected from a human subject while performing mind wandering task (Grandchamp et al., 2014). In the first phase, the network encodes N number of frequency components of an EEG signal collected from one channel of EEG recording during the mentioned experiment. There is a constraint imposed on the magnitude of the lateral power coupling weights as shown in Figure 14. This constraint allows only oscillators with nearby frequency components to interact with each other as initially, the natural frequencies of these oscillators are sorted at an increasing order after sampling from a uniform probability distribution ranging between 0 and 5 Hertz.
In the second phase, the network tries to reconstruct the EEG signals collected simultaneously from the other channels during the same experiment. A set of EEG signals from 5 other channels are reconstructed by the network. The reconstruction accuracy of the modeled EEG signals depends on the number of oscillators in the oscillatory layer and the number of desired signals to be reconstructed. Using complex weights for the feed-forward connections is a key factor as it helps to learn the amplitude, as well as phase of the Fourier decomposition of the EEG signal to be reconstructed, provided the frequency components are already learnt in the encoding phase.

DISCUSSION
The unique contribution of the proposed network of complex neural oscillators is the notion of power coupling by which it becomes possible to achieve a stable normalized phase relationship between two oscillators with arbitrary natural frequencies. With such a feature as the backbone, it is possible to construct a network of oscillators that can learn to reconstruct multiple time series signals. Another positive feature of the proposed model is the biological feasibility of the learning mechanisms. The lateral connections in the oscillatory reservoir FIGURE 14 | (A) The reconstructed signal by the network of 100 oscillators after 30 epochs, the teaching signal, and the respective frequency spectrum in the first phase of training. In this phase of training, the D (t) signal is an EEG signal of duration 10 s, and the exactly similar EEG signal is presented for consecutive 30 epochs. The learning parameters are: A ij = 10 −1 for |j − i| < 5 and A ij = 0 for |j − i| ≥ 5 η ω = 0.1, η α = 10 −4 , τ W = 10 4 , ε = 0.5, N = 100, dt = 10 −3 sec n epoch = 30. (B) After the second phase of training, the network reconstructed and the original EEG signals collected from 5 channels during the same experiment from which D (t) for the first phase of training was collected and the corresponding frequency spectrum. Learning parameters: η K = 3 × 10 −5 , η ζ = 10 −6 , no. of epochs for the batch mode of learning = 1000. are trained by a Hebb-like rule, while the forward connections are trained by a modified delta rule adapted to the complex domain. Experimental observations reveal that local interaction occurs at the single neuron level through single synapse whereas longrange interaction occurs between various cortical areas mediated by multiple synapses (Supplementary Figure 6). At single-cell level the synaptic property can define the relative phase difference between tonically firing presynaptic and postsynaptic neurons. Neuronal population-level interaction between various cortical areas can facilitate cross-frequency coupling ensuring a certain phase relationship (Hyafil et al., 2015) between the LFP activity of different frequencies of the respective cortical areas.

Other Complex Valued Network Models
There is a vast literature pertaining to complex-valued neural network models, which are produced by extending their realvalued counterparts to the complex domain. Complex feedforward networks trained by the complex backpropagation algorithm, complex Hopfield network, complex self-organizing maps are examples of such models. Complex formalism plays a key role in linear systems theory in which both systems and signals can be represented as complex functions. As a special case of this, in electrical circuit theory, oscillations are modeled as complex numbers called phasors, that represent the phase of the oscillations. Therefore, there is a longstanding association between complex numbers and oscillations. Nevertheless, complex-valued neural network literature did not seem to have exploited this association. For the most part, complex-valued neural networks operate like 2n dimensional versions of their n-dimensional real-valued counterparts. The proposed model attempts to take an important step in filling this lacuna. By describing oscillations in the complex plane, and invoking the power coupling principle, it is able to elegantly overcome some of the difficulties involved in harnessing the dynamics of multioscillator networks with real coupling.
Network models of complex-valued oscillators have been described before. For example, subsequent studies by Izhikevich (1996, 2000) have shown that weakly connected network of complex oscillators with identical natural frequencies and self-adjoint matrix of complex coupling coefficient can achieve any arbitrary phase difference between two oscillators. It has also been shown that such network can store and retrieve at least one memory pattern in terms of the phase difference of the oscillators by enabling the Hebbian-like learning rule of the complex coupling coefficient as given in eq. 5. Recent study by Kim and Large (2021) has extended the steady-state stability analysis of such a pair of oscillators to the special case of frequency detuning (ω 1 = ω 2 ) for a range of parameter space. In the same study (Kim and Large, 2021), the authors have proposed a generalized framework for a pair of coupled multifrequency complex oscillators (Hoppensteadt and Izhikevich, 1997) with a Hebbian like learning rule for the complex coupling coefficient. Unlike the proposed power coupling strategy, the coupling is defined for oscillators with natural frequencies near the resonant condition: mω 1 = nω 2 , where m and n are integer numbers. Similar steady-state stability analysis of such a pair of complex oscillators (for n = 2, m = 1 and higher order) shows that the angle of the complex coupling coefficient (θ ij ) learns m∅ 1 − n∅ 2 for 12 = mω 1 − nω 2 = 0, and increases/decreases linearly or rotates for 12 = 0 at steady state. However, neither of the study extended the result to the case of interacting oscillators with very different natural frequencies (|mω 1 − nω 2 | > , where is a small positive number), as we have done using power coupling. Table 1 summarizes the comparison between four different coupling strategies between a pair of complex Hopf oscillators. Comparison is made principally w.r.t two dynamical properties: entrainment and synchronization and Hebbian like plasticity of the coupling coefficient and the coupling term influencing the oscillator dynamics. Frequency entrainment is the dynamical phenomena when the natural frequency of two coupled autonomous oscillators converges to an intermediate frequency (Pikovsky et al., 2001). Whereas the synchronization has been defined here in a more general sense: any coupled oscillators can exhibit synchronization if they maintain any of the phase relationships ( ω j ) constant at steady-state with or without being entrained. A network of complex-valued oscillators that can store patterns as oscillatory states was described in Chakravarthy and Ghosh (1996). The model also proposed a complex form of Hebb's rule, similar to the one used in the current study. However, the model of Chakravarthy and Ghosh (1996) was limited in its capability by the fact that all the oscillators in the model have a common frequency. A representation of synaptic strength using a complex weight, instead of the usual real weight, has an added advantage in representing the temporal relationships underlying neural dynamics. We have recently shown (Chakravarthy, 2020) that the imaginary part of the complex weight captures the temporal asymmetry between the activities of pre-and post-synaptic neurons in a manner akin to the weight kernel of Spike Time Dependent Plasticity (STDP) mechanism (Bi and Poo, 1998).
It may be said that the network of Section "A Generative Network Which Is Capable of Modeling EEG Signals" achieves a discrete approximation of the continuous spectra of the signals that are reconstructed by the network. The greater the number of oscillators in the network, the tighter is the approximation. In this approximation, the trainable real feed-forward weights (α i ) of the network learn the amplitude spectrum of D (t). The angle of the power coupling weight learns normalized phase differences,

Analogous Trainable Central Pattern Generator Model and Comparison
It has previously been shown in Righetti et al. (2005) that, when a network of supercritical Hopf oscillators driven by the reconstruction error signal, individual oscillators are tuned to the nearest frequency components in the target time series. In the network (net-1: schematic is described in Figure 2, and the network dynamics and the learning dynamics is given in eqs. 4 to 8 of Righetti et al., 2005) one shortcoming is that there is no way to recreate the teaching time-series signal (P teach (t)) back without the error feedback loop as with the learnt α i s and ω i s, when the network is reset (i.e., r i s Here synchronization is used in the general sense, i.e., two oscillations are said to be synchronized if either any of the following terms are constant: . RC: Real coupling strategy, CC: complex coupling strategy (Hoppensteadt and Izhikevich, 1996), MFCC: multi-frequency complex coupling strategy (Kim and Large, 2021), CPC: complex power coupling strategy. : possible and/or defined, : possible and/or defined under constrained parameter space, : not possible and/or undefined. Subscripts i and j represent presynaptic and postsynaptic oscillators respectively. and ∅ i s are reset) there is no external forcing to drive the phase offset of each oscillator to the desired phase offsets ϕ i s . To overcome this issue, the authors have proposed a phase signal of frequency same as i th oscillator but same phase as the oscillator with the lowest frequency as a signal communicated by the oscillator with the lowest frequency to the i th oscillator, which makes those two oscillators maintain a certain normalized phase difference ∅ i ω i − ∅ 0 ω 0 . However, to learn/store the phase relationship between the lowest frequency component and the i th frequency component in the teaching signal a "phase variable" considered as a third variable of the oscillator dynamics has been defined which learns the difference between phase signal from the oscillator with the lowest frequency and its own phase. The main issue with this network is that the network can only learn Fourier decomposition (frequency, magnitude and relative phase offset among the constituting frequency components) of a periodic signal which has frequency components of the form ω i = nω 0 , ω 0 being the fundamental frequency of the periodic signal as n ∈ N as the phase signal R i from the 0 th oscillator is interacting with the x dynamics of the 0 th oscillator through real coupling constant. Furthermore, there is an ambiguity in defining the lowest frequency component or fundamental frequency of the periodic teaching signal as the dynamics do not itself identify the oscillator, which learns the lowest frequency component in the teaching signal.
The above arrangement of representing phase difference using an explicit phase parameter located in the i th oscillator has another drawback of Righetti et al. (2005). Phase difference between two oscillators is a property that arises out of their interactions. Ideally, it must be represented as a property of the interaction strength, which in our case is the complex weight parameter. In Righetti et al. (2005), since phase difference information between two oscillators has to be encoded within one of the oscillators, in a rather awkward fashion, each oscillator is restricted to encode the phase difference only with one "standard" oscillator -the oscillator corresponding to the fundamental frequency. On the contrary, in our proposed model, (normalized) phase difference arises jointly out of the power coupling and the corresponding complex weight, there is no need for the abovementioned restrictions. The oscillators can independently control their (normalized) phase differences with other oscillators via their mutual coupling parameters.
To overcome these issues, the proposed network as described in section "A Generative Network Which Is Capable of Modeling EEG Signals" is a similar network where the reservoir of complex supercritical Hopf oscillators is coupled through power coupling to preserve the phase relationship among the constituting frequency components in the teaching signal. When second network architecture proposed by Righetti et al. (schematic of which is described in Figure 3 and the dynamics is given in eqs. 13-20 of Righetti et al. (2005) tries to learn Fourier decomposition of a teaching signal the natural frequencies of the oscillators and the real feed-forward weights learns the magnitude spectrum of the teaching signal. During learning the error signal, F (t) drives the phase offset (∅ i −ω i t) of each of the oscillators to the desired phase offset (ϕ i ) can be retrieved from the phase spectrum of the teaching signal. Thus, the same network deliberated by eq. 20 with the oscillators coupled through power coupling with fixed A ij ( µ) ensures θ ij to learn ϕ i ω j −ϕ j ω i .
One issue while retrieving the stored oscillatory pattern is the dependence on the initial state of the oscillators to achieve the desired solution (see section "Hebbian Learning for the N-Oscillator System With Power Coupling"). We observed that when ∅ i (0) s are sampled from space ∧, subspace of ∧ = {x ∈ R N |0 ≤ x ≤ 2π}, the network attains the desired solution at a steady-state. However, there are other plausible solutions or spurious states the network can acquire. The desired solution of σ ij * (=0) and the spurious solutions are independent of the network parameters like A ij (for A ij = A) and θ ij but has a dependence on the natural frequency of the oscillators. However, the solutions of ψ ij are dependent on both network parameters (θ ij and ω i , assuming A ij = A). The dependency of the boundary of space ? or the boundary of the other subspaces (initializing ∅ i by sampling which leads the system to spurious states) on the natural frequency of the oscillators is yet to be understood. A brief comparison between aCPG model of Righetti et al. (2005) and the proposed aCPG model is summarized in the following Table 2.

Large-Scale Network Models
In the present study, we use the oscillatory neural network model to reconstruct six EEG time series simultaneously. However, the proposed learning mechanisms and network architecture can be easily scaled up. It is possible to model a much larger number of EEG channels (say, 128 or 512) simply by increasing the size of the oscillatory reservoir and the number of output neurons. In the current study, the connections among the oscillators are constrained to ensure that only oscillators with nearby natural frequencies are connected. In a large-scale model, it is conceivable to associate the individual neural oscillators with anatomical locations in the brain, and constrain their coupling based on structural connectivity information from structural imaging tools like Diffusion Tensor Imaging (DTI) (Damoiseaux and Greicius, 2009). Through such developments, the present model can evolve into a class large-scale models of brain dynamics similar to The Virtual Brain (TVB) model, with some relative advantages.
The TVB model principally uses neuronal mass models like the two-variable nonlinear oscillatory models including Fitzhugh-Nagumo model, Winson-Cowan model, Wong-Wang model, Brunel-Wang model, Jensen-Rit model, Stefanescu-Jirsa model (Wilson and Cowan, 1972;Jansen and Rit, 1995;Brunel and Wang, 2003;Wong and Wang, 2006;Stefanescu and Jirsa, 2008), which can exhibit excitable dynamics and limit cycle oscillations. In these models, the oscillators are coupled by only one variable, which is analogous to "real-valued" coupling. They are trained by global optimization algorithms that adjust the coupling strengths so that the output of the oscillatory network matches the recorded brain dynamics. Furthermore, most of them do not have an explicit natural frequency as an independent parameter which is trainable depending on the external input signal. Powanwe and Longtin (2019) developed a large-scale cortical model using a network of stochastic Wilson-Cowan model with 'smooth nonlinearity' to model gamma-band LFP activity. The LFP activity of excitatory and inhibitory population of neurons is modeled using sinusoidal variations (V x (t) = Z x (t) cos (ω 0 t + ∅ x (t))). The three characteristic features, envelop (Z x (t)), constant mean frequency (ω 0 ) and phase (∅ x (t)) respectively are dependent on network parameters. ω 0 is dependent on excitatory synaptic strength, tunable through Hebbian plasticity. Other large-scale modeling studies on the same lines (Mejias et al., 2016;Joglekar et al., 2018) focus on investigating the principal network parameters behind signal propagation or intercommunication between various cortical areas located at the various levels of hierarchy. Although these models can capture the large-scale network entrainment by the external input signal through Hebb-like plasticity of the excitatory synaptic connection, they fail to provide simplified network dynamics essential to learn three constitutive features (frequency, magnitude, and phase offset) of any input signal.

Attractor Dynamics
The proposed network in section "Power Coupling Among N Hopf Oscillators" can reproduce a given time series signal as a linear summation of the oscillations with the same frequency components as the natural frequency of the oscillators in the reservoir and the relative phase relationships encoded in terms of the angle of the power coupling. Given the similarities between conventional attractor neural networks such as the Hopfield network, storing multiple "patterns" can be regarded as multiple time series signals with identical frequency components but non-identical relative phase relationships among them. As the relative phases among the oscillations is encoded in terms of the angle of the power coupling coefficient, the p th pattern can be regarded as: x p = e iθ p where, θ p is a N-dimensional vector, N being the no of oscillators.  Righetti et al. (2005) with our proposed aCPG model.