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, 2003, 2004; Ghosh-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 (Wang and 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₁ and x₂) 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.

A Network of Complex-Valued Oscillators With Complex Coupling and Power Coupling

The canonical model of Hopf oscillator without any external input is described as follows in complex state variable representation $(statevariable:z=r{e}^{i\mathrm{\varnothing }};;i=\sqrt{-1})$, Cartesian coordinate representation (state variables: x, y) and polar coordinate representation (state variables: r, ∅) respectively:

Complex state variable representation:

$\stackrel{.}{z}=z\left(\mathrm{\alpha }+i\mathrm{\omega }+{\mathrm{\beta }}_1{\left|z\right|}^2+\frac{\mathrm{ϵ}{\mathrm{\beta }}_2{\left|z\right|}^5}{1-\mathrm{ϵ}{\left|z\right|}^2}\right)$

Cartesian coordinate representation:

$\stackrel{.}{x}=x\left(\mathrm{\alpha }+{\mathrm{\beta }}_1{r}^2+\frac{\mathrm{ϵ}{\mathrm{\beta }}_2{r}^5}{1-\mathrm{ϵ}{r}^2}\right)-\mathrm{\omega }y$

$\stackrel{.}{y}=y\left(\mathrm{\alpha }+{\mathrm{\beta }}_1{r}^2+\frac{\mathrm{ϵ}{\mathrm{\beta }}_2{r}^5}{1-\mathrm{ϵ}{r}^2}\right)+\mathrm{\omega }x$

Polar coordinate representation:

$\stackrel{.}{r}=\mathrm{\alpha }r+{\mathrm{\beta }}_1{r}^3+\frac{\mathrm{ϵ}{\mathrm{\beta }}_2{r}^5}{1-\mathrm{ϵ}{r}^2}$

$\stackrel{.}{\mathrm{\varnothing }}=\mathrm{\omega }$

Depending on the parameter values (α, β₁, β₂), the autonomous behavior of the oscillator principally falls into four regimes: critical Hopf regime (α = 0, β₁ < 0, β₂ = 0), supercritical Hopf regime (α > 0, β₁ < 0, β₂ = 0), supercritical double limit cycle regime (α < 0, β₁ > 0, β₂ < 0, localmaxima > 0) and subcritical double limit cycle regime (α < 0, β₁ > 0, β₂ < 0, localmaxima < 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ω₀t). The supercritical Hopf regime has an unstable fixed point at the origin and a stable limit cycle at $r=\sqrt{\frac{\mathrm{\alpha }}{\left|{\mathrm{\beta }}_1\right|}}$. 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 (ω − ω₀ < ϵ). Throughout the rest of the paper, we will be using supercritical Hopf regime (with α = μ, β₁ = 1, β₂ = 0), which can be defined as follows:

Complex state variable representation:

$\stackrel{.}{z}=z\left(\mathrm{\mu }+i\mathrm{\omega }-{\left|z\right|}^2\right)$

Cartesian coordinate representation:

$\stackrel{.}{x}=x\left(\mathrm{\mu }-r^2\right)-\mathrm{\omega }y$

$\stackrel{.}{y}=y\left(\mathrm{\mu }-r^2\right)+\mathrm{\omega }x$

Polar coordinate representation:

$\stackrel{.}{r}=\mathrm{\mu }r-r^3$

$\stackrel{.}{\mathrm{\varnothing }}=\mathrm{\omega }$

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).

$\stackrel{.}{z}=z\left(\mathrm{\mu }+i\mathrm{\omega }-{\left|z\right|}^2\right)+\mathrm{\epsilon }{I}_{ext}\left(t\right)(\mathrm{1a})$

$\stackrel{.}{\mathrm{\omega }}=-\mathrm{\epsilon }{I}_{ext}\left(t\right)\sin \mathrm{\varnothing }(\mathrm{1b})$

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 $\sqrt{\mathrm{\mu }}$ for ε = 0. If I_ext(t) = I₀sin⁡(ω₀t + φ) (I₀, ω₀, and φ are the magnitude, frequency and phase offset of the sinusoidal input signal), the natural frequency, ω, adapts to the frequency ω₀.

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:

$\begin{array}{c}\hfill \stackrel{.}{z}=z\left(\mathrm{\mu }+i\mathrm{\omega }-{\left|z\right|}^2\right)+\mathrm{\epsilon }{I}_{ext}\left(t\right)\\ \hfill I_{ext}\left(t\right)=I_0{e}^{i\left({\mathrm{\omega }}_{0}t+\mathrm{\phi }\right)}\end{array}(\mathrm{2a})$

$\stackrel{.}{\mathrm{\omega }}=-\mathrm{\epsilon }\left(real\left(I_{ext}\left(t\right)\right)\sin \mathrm{\varnothing }-img\left(I_{ext}\left(t\right)\right)\cos \mathrm{\varnothing }\right)(\mathrm{2b})$

$\stackrel{.}{\mathrm{\omega }}=-\mathrm{\epsilon }{I}_0\sin \left(\mathrm{\varnothing }-{\mathrm{\omega }}_{0}t-\mathrm{\phi }\right)(\mathrm{2c})$

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:

$\begin{array}{c}\stackrel{.}{x}=\left(\mathrm{\mu }-r^2\right)x-\mathrm{\omega }y+\mathrm{\epsilon }{I}_0\cos ({\mathrm{\omega }}_{0}t+\mathrm{\phi })\hfill \\ \stackrel{.}{y}=\left(\mathrm{\mu }-r^2\right)y+\mathrm{\omega }x+\mathrm{\epsilon }{I}_0\sin ({\mathrm{\omega }}_{0}t+\mathrm{\phi })\hfill \end{array}(\mathrm{2a1})$

$\begin{array}{c}\stackrel{.}{r}=\left(\mathrm{\mu }-r^2\right)r+\mathrm{\epsilon }{I}_0\cos \left({\mathrm{\omega }}_{0}t+\mathrm{\phi }-\mathrm{\varnothing }\right)\hfill \\ \stackrel{.}{\mathrm{\varnothing }}=\mathrm{\omega }+\frac{\mathrm{\epsilon }{I}_0}{r}\sin \left({\mathrm{\omega }}_{0}t+\mathrm{\phi }-\mathrm{\varnothing }\right)\hfill \end{array}(\mathrm{2a2})$

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 $\frac{\mathrm{\epsilon }{I}_0}{r}\sin \left({\mathrm{\omega }}_{0}t-\mathrm{\varnothing }\right)$. 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 ($\overrightarrow{e_{\mathrm{\varnothing }}}$), $\overrightarrow{P_{\mathrm{\varnothing }}}$ is (to understand the vector notations refer to Figure 1):

$\begin{array}{c}r\overrightarrow{P_{\mathrm{\varnothing }}}=real\left(\mathrm{\epsilon }{I}_{ext}\left(t\right)\right)\sin \mathrm{\varnothing }-img\left(\mathrm{\epsilon }{I}_{ext}\left(t\right)\right)\cos \mathrm{\varnothing }\hfill \\ =\mathrm{\epsilon }{I}_0\cos \left({\mathrm{\omega }}_{0}t+\mathrm{\phi }\right)\sin \mathrm{\varnothing }-\mathrm{\epsilon }{I}_0\sin \left({\mathrm{\omega }}_{0}t+\mathrm{\phi }\right)\cos \mathrm{\varnothing }\hfill \\ \overrightarrow{P_{\mathrm{\varnothing }}}={\displaystyle \frac{\mathrm{\epsilon }{I}_0}{r}}\sin \left(\mathrm{\varnothing }-{\mathrm{\omega }}_{0}t-\mathrm{\phi }\right)\hfill \end{array}$

FIGURE 1

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).

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.

FIGURE 2

Figure 2. Equations 2a-c is simulated for $\mathrm{\mu }=1,\mathrm{\omega }\left(0\right)=40,{\mathrm{\omega }}_0=30,\mathrm{\epsilon }=0.9,I_0=1,\mathrm{\phi }=\frac{\mathrm{\pi }}{4},dt=0.001sec$ for 1000 s. (A) Depicts the variation of y(t) w.r.t img(I_ext(t)) at various time instants. (B) Depicts the variation of x(t) w.r.t real(I_ext(t)) at various time instants. (C) It can be observed that the natural frequency of the Hopf oscillator adapts to the frequency of the input signal.

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)

$\stackrel{.}{z_1}=z_1\left(\mathrm{\mu }+i{\mathrm{\omega }}_1-{\left|z_1\right|}^2\right)+W_{12}real\left(z_2\right)(\mathrm{3a})$

$\stackrel{.}{z_2}=z_2\left(\mathrm{\mu }+i{\mathrm{\omega }}_2-{\left|z_2\right|}^2\right)+W_{21}real\left(z_1\right)(\mathrm{3b})$

where W₁₂ is the real coupling coefficient from 2nd oscillator to 1st one, and W₂₁ is the real coupling coefficient from 1st to 2nd. In the stated scenario, both of the oscillators have identical natural frequencies, ω₁ = ω₂. 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 coefficient. (for proof refer to the Supplementary Appendix 2).

$\stackrel{.}{z_1}=z_1\left(\mathrm{\mu }+i\mathrm{\omega }-{\left|z_1\right|}^2\right)+W{z}_2(\mathrm{4a})$

$\stackrel{.}{z_2}=z_2\left(\mathrm{\mu }+i\mathrm{\omega }-{\left|z_2\right|}^2\right)+W^{*}{z}_1(\mathrm{4b})$

where z₁ = r₁e^i∅₁, z₂ = r₂e^i∅₂ and W = Ae^iθ, W^∗=Ae^−iθ to be the coupling coefficient (A and θ being the magnitude and the angle of complex coupling coefficient), in polar coordinate system representation:

$\begin{array}{c}\stackrel{.}{r_1}=\left(\mathrm{\mu }-r_1{}^2\right)r_1+A{r}_2\cos \left({\mathrm{\varnothing }}_2-{\mathrm{\varnothing }}_1+\mathrm{\theta }\right)\hfill \\ \stackrel{.}{{\mathrm{\varnothing }}_1}=\mathrm{\omega }+A\frac{r_2}{r_1}\sin ({\mathrm{\varnothing }}_2-{\mathrm{\varnothing }}_1+\mathrm{\theta })\hfill \end{array}(\mathrm{4a1})$

$\begin{array}{c}\stackrel{.}{r_2}=\left(\mathrm{\mu }-r_2{}^2\right)r_2+A{r}_1\cos \left({\mathrm{\varnothing }}_1-{\mathrm{\varnothing }}_2-\mathrm{\theta }\right)\hfill \\ \stackrel{.}{{\mathrm{\varnothing }}_2}=\mathrm{\omega }+A\frac{r_1}{r_2}\sin ({\mathrm{\varnothing }}_1-{\mathrm{\varnothing }}_2-\mathrm{\theta })\hfill \end{array}(\mathrm{4b1})$

At steady-state ∅₁−∅₂ approaches any of the solutions 2nπ + θ depending on the initial conditions [∅₁(0) and ∅₂(0)], whereas the magnitude of the complex coupling coefficient determines the rate at which the phase-locking occurs. ψ_ss (where, ψ = ∅₁−∅₂ and ψ_ss is the steady state value of ψ) attains solution 2nπ + θ for the following initial condition 2nπ−(−π + θ) < ∅₁(0)−∅₂(0) < 2nπ−(π + θ).

Simulated Result

To verify the above result, we have simulated eqs. 4a1 and 4b1 numerically by Euler integration method on MATLAB platform with △t = 1msec. 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.

FIGURE 3

Figure 3. Two Hopf oscillators, dynamics of which is defined in the eqs. 4a1 and 4b1 simulated with ω = 5, μ = 1, A = 0.5, $\mathrm{\theta }=\frac{\mathrm{\pi }}{4}$ for various initial conditions (−π < ∅₁(0)−∅₂(0)−θ≤3π) depicting that ${\mathrm{\psi }}_{ss}\left(=\mathrm{\varnothing }_1{}_{ss}-\mathrm{\varnothing }_2{}_{ss}\right)$ can reach any of the following solutions 2nπ + θ depending on the initial condition.

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.

To train the complex coupling coefficient, a Hebbian like learning rule can be used as follows:

${\mathrm{\tau }}_W\stackrel{.}{W}=-W+z_{1}z_2{}^{*}(5)$

The polar coordinate representation:

${\mathrm{\tau }}_W\stackrel{.}{A}=-A+r_1{r}_2\cos \left({\mathrm{\varnothing }}_1-{\mathrm{\varnothing }}_2-\mathrm{\theta }\right)(\mathrm{5a})$

${\mathrm{\tau }}_W\stackrel{.}{\mathrm{\theta }}=\frac{r_1{r}_2}{A}\sin ({\mathrm{\varnothing }}_1-{\mathrm{\varnothing }}_2-\mathrm{\theta })(\mathrm{5b})$

We have limited the scope of our study only to θ dynamics by assuming $\stackrel{.}{A}=0$. Assuming τ_W is the time constant for the specified learning dynamics. A similar adaptation of the Hebb’s rule was used earlier to train connections in a complex version of the Hopfield network (Hopfield, 1982). When the dynamics of two Hopf oscillators is defined by eq. 4 (with no external input), and the complex coupling coefficient is trained according to eq. 5 with very small value of A(|W|), θ(angle(W)) learns ∅₁(0) − ∅₂(0). From eq. 4, it can be interpreted that ∅₁(t) ≈ ωt + ∅₁(0) and ∅₂(t) ≈ ωt + ∅₂(0)(Assuming A≪1), similarly from eq. 5, we can see that θ or angle(W) learns $angle\left({\text{z}}_1\text{z}_2{}^{*}\right)$ or ∅₁(t) − ∅₂(t) ≈ ∅₁(0) − ∅₂(0). It is shown in the Supplementary Appendix 3 that the nonzero symmetric steady state solution (r_1ss = r_2ss ≠ 0, θ − ψ_ss) is $r_{1ss}=\sqrt{\mathrm{\mu }+A}$, ψ = 2nπ + θ, which is a stable node.

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 (φ₁, φ₂) will drive those oscillators to oscillate at the same phase as the external sinusoidal forcing (∅₁ ≈ ωt + φ₁ and ∅₂ ≈ ωt + φ₂), provided very low magnitude of the complex coupling coefficient (A≪1).

$\begin{array}{c}\stackrel{.}{z_1}=z_1\left(\mathrm{\mu }+i\mathrm{\omega }-{\left|z_1\right|}^2\right)+A{e}^{i\mathrm{\theta }}{z}_2+I_{01}{e}^{i(\mathrm{\omega }t+{\mathrm{\phi }}_1)}\hfill \\ \stackrel{.}{z_2}=z_2\left(\mathrm{\mu }+i\mathrm{\omega }-{\left|z_2\right|}^2\right)+A{e}^{-i\mathrm{\theta }}{z}_1+I_{02}{e}^{i(\mathrm{\omega }t+{\mathrm{\phi }}_2)}\hfill \end{array}(6)$

Under the influence of external input, the angles of z₁ and z₂ (∅₁ and ∅₂ respectively) tend to ωt + φ₁ and ωt + φ₂, respectively. When the complex coupling coefficient is trained according to eq. 5, the angle(W) or θ learns $angle\left({\text{z}}_1\text{z}_2{}^{*}\right)$ or ∅₁(t)−∅₂(t) ≈ φ₁−φ₂. During training, A must be kept low (A≪1) so that z₁ and z₂ 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 $({\mathrm{\phi }}_1-\mathrm{\phi }{}_2)$.

FIGURE 4

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 φ₁−φ₂. I₀₁ = I₀₂ = 0.3, ${\mathrm{\phi }}_1=\frac{\mathrm{\pi }}{4},{\mathrm{\phi }}_2=\frac{\mathrm{\pi }}{6},\mathrm{\omega }=5,A={10}^{-5},\mathrm{\eta }=\frac{1}{dt},dt=0.001$.

Novel Contributions: Power Coupling

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 natural frequencies, depending on the coupling strength.

${\stackrel{.}{\mathrm{\varnothing }}}_1={\mathrm{\omega }}_1+K_1\sin ({\mathrm{\varnothing }}_2-{\mathrm{\varnothing }}_1)$

${\stackrel{.}{\mathrm{\varnothing }}}_2={\mathrm{\omega }}_2+K_2\sin ({\mathrm{\varnothing }}_1-{\mathrm{\varnothing }}_2)$

Let, K₁ and K₂ be the real coupling strengths in the forward and backward directions, respectively. The phase difference between the two oscillators (∅ = ∅₁ − ∅₂) reaches a steady-state value if K₁ + K₂ > |ω₁ − ω₂| which is satisfied by the relationship sin⁡∅^∗=$\frac{{\mathrm{\omega }}_1-{\mathrm{\omega }}_2}{K_1+K_2}$ where ∅^∗ be the steady-state phase difference. Notably, both of the oscillators reach a periodic solution with frequency ω^∗=$\frac{K_1{\mathrm{\omega }}_2+K_2{\mathrm{\omega }}_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).

${\stackrel{.}{\mathrm{\varnothing }}}_i={\mathrm{\omega }}_i+\frac{K}{N}\sum _{j=1}^N\sin \left({\mathrm{\varnothing }}_j-{\mathrm{\varnothing }}_i\right)(\mathrm{7a})$

$\mathrm{\hspace{1em}}\mathit{\hspace{1em}\hspace{1em}}r{e}^{i\mathrm{\theta }}=\frac{1}{N}\sum _{j=1}^N{e}^{i{\mathrm{\varnothing }}_j}$

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 (∅₁ − ∅₂), although when the coupled oscillators are entrained in the m:n ratio, the phase difference has been defined as m∅₁ − n∅₂ (where, m, n are integers). Therefore, one may define a quantity called “normalized phase difference” as follows:

${\mathrm{\psi }}_{12}=-{\mathrm{\psi }}_{21}=\frac{{\mathrm{\varnothing }}_1}{{\mathrm{\omega }}_1}-\frac{{\mathrm{\varnothing }}_2}{{\mathrm{\omega }}_2}(\mathrm{7b})$

Let ψ₁₂ and ψ₂₁ 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:

$\stackrel{.}{z_1}=z_1\left(\mathrm{\mu }+i{\mathrm{\omega }}_1-{\left|z_1\right|}^2\right)+A_{12}{e}^{i\frac{{\mathrm{\theta }}_{12}}{{\mathrm{\omega }}_2}}z_2{}^{\frac{{\mathrm{\omega }}_1}{{\mathrm{\omega }}_2}}(\mathrm{8a})$

$\stackrel{.}{z_2}=z_2\left(\mathrm{\mu }+i{\mathrm{\omega }}_2-{\left|z_2\right|}^2\right)+A_{21}{e}^{i\frac{{\mathrm{\theta }}_{21}}{{\mathrm{\omega }}_1}}z_1{}^{\frac{{\mathrm{\omega }}_2}{{\mathrm{\omega }}_1}}(\mathrm{8b})$

where, $A_{12}{e}^{i\frac{{\mathrm{\theta }}_{12}}{{\mathrm{\omega }}_2}}=W_{12}$ is the weight of power coupling from 2nd oscillator to the 1st oscillator and $A_{21}{e}^{i\frac{{\mathrm{\theta }}_{21}}{{\mathrm{\omega }}_1}}=W_{21}$ is the weight of power coupling from 1st oscillator to the 2nd oscillator.

The polar coordinate representation:

$\stackrel{.}{r_1}=\left(\mathrm{\mu }-r_1{}^2\right)r_1+A_{12}r_2{}^{\frac{{\mathrm{\omega }}_1}{{\mathrm{\omega }}_2}}\cos {\mathrm{\omega }}_1\left({\displaystyle \frac{{\mathrm{\varnothing }}_2}{{\mathrm{\omega }}_2}}-{\displaystyle \frac{{\mathrm{\varnothing }}_1}{{\mathrm{\omega }}_1}}+{\displaystyle \frac{{\mathrm{\theta }}_{12}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}}\right)\hspace{1em}\text{(8a1)}$

$\stackrel{.}{{\mathrm{\varnothing }}_1}={\mathrm{\omega }}_1+A_{12}{\displaystyle \frac{r_2{}^{\frac{{\mathrm{\omega }}_1}{{\mathrm{\omega }}_2}}}{r_1}}\sin {\mathrm{\omega }}_1\left({\displaystyle \frac{{\mathrm{\varnothing }}_2}{{\mathrm{\omega }}_2}}-{\displaystyle \frac{{\mathrm{\varnothing }}_1}{{\mathrm{\omega }}_1}}+{\displaystyle \frac{{\mathrm{\theta }}_{12}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}}\right)\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.25em}\text{(8a2)}$

$\stackrel{.}{r_2}=\left(\mathrm{\mu }-r_2{}^2\right)r_2+A_{21}r_1{}^{\frac{{\mathrm{\omega }}_2}{{\mathrm{\omega }}_1}}\cos {\mathrm{\omega }}_2\left({\displaystyle \frac{{\mathrm{\varnothing }}_1}{{\mathrm{\omega }}_1}}-{\displaystyle \frac{{\mathrm{\varnothing }}_2}{{\mathrm{\omega }}_2}}+{\displaystyle \frac{{\mathrm{\theta }}_{21}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}}\right)\text{(8b1)}$

$\stackrel{.}{{\mathrm{\varnothing }}_2}={\mathrm{\omega }}_2+A_{21}{\displaystyle \frac{r_1{}^{\frac{{\mathrm{\omega }}_2}{{\mathrm{\omega }}_1}}}{r_2}}\sin {\mathrm{\omega }}_2\left({\displaystyle \frac{{\mathrm{\varnothing }}_1}{{\mathrm{\omega }}_1}}-{\displaystyle \frac{{\mathrm{\varnothing }}_2}{{\mathrm{\omega }}_2}}+{\displaystyle \frac{{\mathrm{\theta }}_{21}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}}\right)\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.5em}\hspace{0.25em}\text{(8b2)}$

The schematic of the coupling architecture is elaborated in the following Figure 5.

FIGURE 5

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:$\frac{{\mathrm{\omega }}_j}{{\mathrm{\omega }}_i}$. 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}=A{e}^{i\frac{{\mathrm{\theta }}_{ij}}{{\mathrm{\omega }}_j}}$ to convert the forcing signal from j^th oscillator to i^th oscillator and vice versa.

Based on eqs. 8a2 and 8b2, we understand that two modified Kuramoto oscillators may be coupled using power coupling as follows:

$\begin{array}{c}\stackrel{.}{{\mathrm{\varnothing }}_1}={\mathrm{\omega }}_1+A_{12}\sin {\mathrm{\omega }}_1\left(\frac{{\mathrm{\varnothing }}_2}{{\mathrm{\omega }}_2}-\frac{{\mathrm{\varnothing }}_1}{{\mathrm{\omega }}_1}+\frac{{\mathrm{\theta }}_{12}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}\right)\hfill \\ \stackrel{.}{{\mathrm{\varnothing }}_2}={\mathrm{\omega }}_2+{\text{A}}_{21}\sin {\mathrm{\omega }}_2\left(\frac{{\mathrm{\varnothing }}_1}{{\mathrm{\omega }}_1}-\frac{{\mathrm{\varnothing }}_2}{{\mathrm{\omega }}_2}+\frac{{\mathrm{\theta }}_{21}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}\right)\hfill \end{array}(9)$

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:

$\stackrel{.}{z_i}=z_i\left(\mathrm{\mu }+i{\mathrm{\omega }}_i-{\left|z_i\right|}^2\right)+\sum _{j=1\ni j\ne i}^N{A}_{ij}{e}^{i\frac{{\mathrm{\theta }}_{ij}}{{\mathrm{\omega }}_j}}z_j{}^{\frac{{\mathrm{\omega }}_i}{{\mathrm{\omega }}_j}}(10)$

Let $W_{ij}=A_{ij}{e}^{i\frac{{\mathrm{\theta }}_{ij}}{{\mathrm{\omega }}_j}}$ be the weight of power coupling from the j^th oscillator to i^th oscillator.

The polar coordinate representation is:

$\stackrel{.}{r_i}=\left(\mathrm{\mu }-r_i{}^2\right)r_i+\sum _{j=1\ni j\ne i}^N{A}_{ij}r_j{}^{\frac{{\mathrm{\omega }}_{\text{i}}}{{\mathrm{\omega }}_{\text{j}}}}\cos {\mathrm{\omega }}_{\text{i}}\left(\frac{{\mathrm{\varnothing }}_j}{{\mathrm{\omega }}_j}-\frac{{\mathrm{\varnothing }}_i}{{\mathrm{\omega }}_i}+\frac{{\mathrm{\theta }}_{ij}}{{\mathrm{\omega }}_{\text{i}}{\mathrm{\omega }}_{\text{j}}}\right)(\mathrm{11a})$

$\stackrel{.}{{\mathrm{\varnothing }}_i}={\mathrm{\omega }}_i+\sum _{j=1\ni j\ne i}^N{A}_{ij}\frac{r_j{}^{\frac{{\mathrm{\omega }}_i}{{\mathrm{\omega }}_j}}}{r_i}\sin {\mathrm{\omega }}_i\left(\frac{{\mathrm{\varnothing }}_j}{{\mathrm{\omega }}_j}-\frac{{\mathrm{\varnothing }}_i}{{\mathrm{\omega }}_i}+\frac{{\mathrm{\theta }}_{ij}}{{\mathrm{\omega }}_i{\mathrm{\omega }}_j}\right)(\mathrm{11b})$

Therefore, the dynamics of the i^th oscillator in a network of N number Kuramoto oscillators coupled through power coupling can be represented as:

$\stackrel{.}{{\mathrm{\varnothing }}_i}={\mathrm{\omega }}_i+\sum _{j=1\ni j\ne i}^N{A}_{ij}\sin {\mathrm{\omega }}_i\left(\frac{{\mathrm{\varnothing }}_j}{{\mathrm{\omega }}_j}-\frac{{\mathrm{\varnothing }}_i}{{\mathrm{\omega }}_i}+\frac{{\mathrm{\theta }}_{ij}}{{\mathrm{\omega }}_i{\mathrm{\omega }}_j}\right)(12)$

We may appreciate the difference in the second term in eqs. 7a and 12. $\stackrel{.}{{\mathrm{\psi }}_{ij}}{}^{\prime }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:

$\begin{array}{c}\stackrel{.}{{\mathrm{\psi }}_{12}}=\frac{A_{12}}{{\mathrm{\omega }}_1}\frac{r_2{}^{\frac{{\mathrm{\omega }}_1}{{\mathrm{\omega }}_2}}}{r_1}\sin {\mathrm{\omega }}_1\left(\frac{{\mathrm{\varnothing }}_2}{{\mathrm{\omega }}_2}-\frac{{\mathrm{\varnothing }}_1}{{\mathrm{\omega }}_1}+\frac{{\mathrm{\theta }}_{12}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}\right)\hfill \\ -\frac{A_{21}}{{\mathrm{\omega }}_2}\frac{r_1{}^{\frac{{\mathrm{\omega }}_2}{{\mathrm{\omega }}_1}}}{r_2}\sin {\mathrm{\omega }}_2\left(\frac{{\mathrm{\varnothing }}_1}{{\mathrm{\omega }}_1}-\frac{{\mathrm{\varnothing }}_2}{{\mathrm{\omega }}_2}+\frac{{\mathrm{\theta }}_{21}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}\right)\hfill \end{array}(\mathrm{13a})$

$\begin{array}{c}\stackrel{.}{{\mathrm{\psi }}_{ij}}=\sum _{k=1\ni k\ne i}^N\frac{A_{ik}}{{\mathrm{\omega }}_{\text{k}}}\frac{r_k{}^{\frac{{\mathrm{\omega }}_i}{{\mathrm{\omega }}_k}}}{r_i}\sin {\mathrm{\omega }}_i\left(\frac{{\mathrm{\varnothing }}_k}{{\mathrm{\omega }}_k}-\frac{{\mathrm{\varnothing }}_i}{{\mathrm{\omega }}_i}+\frac{{\mathrm{\theta }}_{ik}}{{\mathrm{\omega }}_i{\mathrm{\omega }}_k}\right)\hfill \\ -\sum _{k=1\ni k\ne j}^N\frac{A_{jk}}{{\mathrm{\omega }}_{\text{k}}}\frac{r_k{}^{\frac{{\mathrm{\omega }}_j}{{\mathrm{\omega }}_k}}}{r_j}\sin {\mathrm{\omega }}_j\left(\frac{{\mathrm{\varnothing }}_k}{{\mathrm{\omega }}_k}-\frac{{\mathrm{\varnothing }}_j}{{\mathrm{\omega }}_j}+\frac{{\mathrm{\theta }}_{jk}}{{\mathrm{\omega }}_j{\mathrm{\omega }}_k}\right)\hfill \end{array}(\mathrm{13b})$

At steady state, ψ₁₂ and ${\mathrm{\psi }}_{ij}\left(=\frac{{\mathrm{\varnothing }}_i}{{\mathrm{\omega }}_i}-\frac{{\mathrm{\varnothing }}_j}{{\mathrm{\omega }}_j}\right)$ 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, ${\mathrm{\psi }}_{ij}^{\prime }$s can achieve the following desirable solution:

$\mathrm{\psi }_{ij}{}^{*}=\frac{{\mathrm{\theta }}_{ij}}{{\mathrm{\omega }}_i{\mathrm{\omega }}_j}(\mathrm{13b.s})$

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 $\mathrm{\psi }_{ij}{}^{*}$ can be achieved only for certain initial values of $\mathrm{\varnothing }{}_i^{\prime }$s.

Two Hopf Oscillators With Power Coupling

Under the following constraints θ₁₂ = θ, θ₂₁ = −θ, A₁₂ = A₂₁ = A, μ = 1, some of the steady-state solutions of eq. 13a where both entities on the right-hand side of eq. 13a is zero, is as follows:

$\mathrm{\psi }_{12}{}_{ss}={\displaystyle \frac{\mathrm{\theta }}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}}\pm {\displaystyle \frac{n_1\mathrm{\pi }}{{\mathrm{\omega }}_1}}={\displaystyle \frac{\mathrm{\theta }}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}}\pm {\displaystyle \frac{n_2\mathrm{\pi }}{{\mathrm{\omega }}_2}}$

following the specification: $\frac{n_1}{n_2}=\frac{{\mathrm{\omega }}_1}{{\mathrm{\omega }}_2}$. It can be noticed that n₁ = n₂ = 0 gives us the desired solution, as mentioned in eq. 13b.s.

Let ${\mathrm{\sigma }}_{12}={\mathrm{\psi }}_{12}-\frac{\mathrm{\theta }}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}=\frac{{\mathrm{\varnothing }}_1}{{\mathrm{\omega }}_1}-\frac{{\mathrm{\varnothing }}_2}{{\mathrm{\omega }}_2}-\frac{\mathrm{\theta }}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}$, so, ${\mathrm{\sigma }}_{ij}=\frac{{\mathrm{\varnothing }}_i}{{\mathrm{\omega }}_i}-\frac{{\mathrm{\varnothing }}_j}{{\mathrm{\omega }}_j}-\frac{{\mathrm{\theta }}_{ij}}{{\mathrm{\omega }}_i\mathrm{\omega }{}_j}$.

From the following Figure 6A it can be verified that ψ₁₂ and σ₁₂ attain $\frac{\mathrm{\theta }}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}=\frac{-1.8968}{5\times 10}\approx -0.038$ and zero solution (the desired solution as mentioned before) respectively for the initialization specified in the figure caption. To check the dependency of $\mathrm{\sigma }_{12}{}_{ss}$ on the initial values of ∅₁ and ∅₂, ∅₁(0) and ∅₂(0) space is discretized at △∅₁(0) = △∅₂(0) = 0.1, in the interval 0 < ∅₁(0), ∅₂(0) ≤ 2π and eqs. 8a1, 8a2, 8b1, 8b2 is simulated for 200 s to ensure the steady-state to be achieved for various combinations of ∅₁(0) and ∅₂(0). From the following Figure 6B it can be observed that $\mathrm{\sigma }_{12}{}_{ss}$ achieves the following solutions: $0,\pm \frac{2\mathrm{\pi }}{{\mathrm{\omega }}_1}$; for n₁ = 0, 1 or $n_2=0,2(=\frac{{\mathrm{\omega }}_2{n}_1}{{\mathrm{\omega }}_1}=2n_1)$ independent of θ.

FIGURE 6

Figure 6. (A) The variation of ψ, $\stackrel{.}{\mathrm{\psi }}$, and $\mathrm{\psi }-\frac{\mathrm{\theta }}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}$ w.r.t time is plotted by simulating eqs. 8a1 and 8b1 for the following set of parameters and initial conditions. ω₁ = 5, ω₂ = 10, A = 0.05, θ = −1.8968, ∅₁(0) = 3.7008, ∅₂(0) = 2.3106. (B) It can be verified that depending on the various combinations of initial conditions (0 < ∅₁(0), ∅₂(0)≤2π), $\mathrm{\sigma }_{12}{}_{ss}$ can reach one of the solutions $\pm \frac{2n_1\mathrm{\pi }}{{\mathrm{\omega }}_1}=\pm \frac{2n_2\mathrm{\pi }}{{\mathrm{\omega }}_2}$. For the following set of parameters, simulated results show that $\mathrm{\sigma }_{12}{}_{ss}$ reaches one of the following three solutions $0,\pm \frac{2\mathrm{\pi }}{{\mathrm{\omega }}_1}or\pm \frac{4\mathrm{\pi }}{{\mathrm{\omega }}_2}or\pm 1.257$. Simulation parameters: ω₁ = 5, ω₂ = 10, A = 0.2, θ = 2.9644.

Three Hopf Oscillators Coupled Through Power Coupling

We now explore the solution space numerically for a three-oscillator system. When three Hopf oscillators are coupled according to the eq. 10 or 11 (N = 3) with ${\mathrm{\theta }}_{ij}^{\prime }$s, chosen as θ_ij = θ_i − θ_j∀ 0 < θ_i, θ_j ≤ 2π (i.e., θ_ij = −θ_ij), at steady state σ_ij achieves one of the possible solutions depending on the initializations (∅_i(0)′s). The $\stackrel{.}{{\mathrm{\varnothing }}_i}$ dynamics of the three oscillators can be expressed in the following form.

$\begin{array}{c}\stackrel{.}{{\mathrm{\varnothing }}_1}={\mathrm{\omega }}_1+A_{12}\frac{r_2{}^{\frac{{\mathrm{\omega }}_1}{{\mathrm{\omega }}_2}}}{r_1}\sin {\mathrm{\omega }}_1\left(\frac{{\mathrm{\varnothing }}_2}{{\mathrm{\omega }}_2}-\frac{{\mathrm{\varnothing }}_1}{{\mathrm{\omega }}_1}+\frac{{\mathrm{\theta }}_{12}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}\right)\hfill \\ +A_{13}\frac{r_3{}^{\frac{{\mathrm{\omega }}_1}{{\mathrm{\omega }}_3}}}{r_1}\sin {\mathrm{\omega }}_1\left(\frac{{\mathrm{\varnothing }}_3}{{\mathrm{\omega }}_3}-\frac{{\mathrm{\varnothing }}_1}{{\mathrm{\omega }}_1}+\frac{{\mathrm{\theta }}_{13}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_3}\right)\hfill \end{array}(\mathrm{14a})$

$\begin{array}{c}\stackrel{.}{{\mathrm{\varnothing }}_2}={\mathrm{\omega }}_2+A_{23}\frac{r_3{}^{\frac{{\mathrm{\omega }}_2}{{\mathrm{\omega }}_3}}}{r_2}\sin {\mathrm{\omega }}_2\left(\frac{{\mathrm{\varnothing }}_3}{{\mathrm{\omega }}_3}-\frac{{\mathrm{\varnothing }}_2}{{\mathrm{\omega }}_2}+\frac{{\mathrm{\theta }}_{23}}{{\mathrm{\omega }}_2{\mathrm{\omega }}_3}\right)\hfill \\ +A_{21}\frac{r_1{}^{\frac{{\mathrm{\omega }}_1}{{\mathrm{\omega }}_3}}}{r_2}\sin {\mathrm{\omega }}_2\left(\frac{{\mathrm{\varnothing }}_1}{{\mathrm{\omega }}_1}-\frac{{\mathrm{\varnothing }}_2}{{\mathrm{\omega }}_2}+\frac{{\mathrm{\theta }}_{21}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}\right)\hfill \end{array}(\mathrm{14b})$

$\begin{array}{c}\stackrel{.}{{\mathrm{\varnothing }}_3}={\mathrm{\omega }}_3+A_{31}\frac{r_1{}^{\frac{{\mathrm{\omega }}_3}{{\mathrm{\omega }}_1}}}{r_3}\sin {\mathrm{\omega }}_3\left(\frac{{\mathrm{\varnothing }}_1}{{\mathrm{\omega }}_1}-\frac{{\mathrm{\varnothing }}_3}{{\mathrm{\omega }}_3}+\frac{{\mathrm{\theta }}_{31}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_3}\right)\hfill \\ +A_{32}\frac{r_2{}^{\frac{{\mathrm{\omega }}_3}{{\mathrm{\omega }}_2}}}{r_3}\sin {\mathrm{\omega }}_1\left(\frac{{\mathrm{\varnothing }}_2}{{\mathrm{\omega }}_2}-\frac{{\mathrm{\varnothing }}_3}{{\mathrm{\omega }}_3}+\frac{{\mathrm{\theta }}_{32}}{{\mathrm{\omega }}_2{\mathrm{\omega }}_3}\right)\hfill \end{array}(\mathrm{14c})$

At steady-state $\stackrel{.}{{\mathrm{\psi }}_{ij}}=0$, where ${\mathrm{\psi }}_{ij}=\frac{{\mathrm{\varnothing }}_i}{{\mathrm{\omega }}_i}-\frac{{\mathrm{\varnothing }}_j}{{\mathrm{\omega }}_j}$ or $\stackrel{.}{{\mathrm{\psi }}_{ij}}=\frac{\stackrel{\_}{{\mathrm{\varnothing }}_i}}{{\mathrm{\omega }}_i}-\frac{\stackrel{\_}{{\mathrm{\varnothing }}_j}}{{\mathrm{\omega }}_j}$: Assuming μ = 1, A_ij = A and e^iθ_ij is a Hermitian matrix (i.e., e^iθ_ij = e^−iθ_ji),

$\begin{array}{c}\stackrel{.}{\mathrm{\psi }_{12}{}_{ss}}=-{\displaystyle \frac{1}{{\mathrm{\omega }}_1}}\sin {\mathrm{\omega }}_1\left(\mathrm{\psi }_{12}{}_{ss}-{\displaystyle \frac{{\mathrm{\theta }}_{12}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}}\right)\hfill \\ +{\displaystyle \frac{1}{{\mathrm{\omega }}_1}}\sin {\mathrm{\omega }}_1\left(\mathrm{\psi }_{31}{}_{ss}-{\displaystyle \frac{{\mathrm{\theta }}_{31}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_3}}\right)+{\displaystyle \frac{1}{{\mathrm{\omega }}_2}}\sin {\mathrm{\omega }}_2\left(\mathrm{\psi }_{23}{}_{ss}-{\displaystyle \frac{{\mathrm{\theta }}_{23}}{{\mathrm{\omega }}_2{\mathrm{\omega }}_3}}\right)\hfill \\ -{\displaystyle \frac{1}{{\mathrm{\omega }}_2}}\sin {\mathrm{\omega }}_2\left(\mathrm{\psi }_{12}{}_{ss}-{\displaystyle \frac{{\mathrm{\theta }}_{12}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}}\right)=0\hfill \end{array}$

$\begin{array}{c}\stackrel{.}{\mathrm{\psi }_{23}{}_{ss}}=-{\displaystyle \frac{1}{{\mathrm{\omega }}_2}}\sin {\mathrm{\omega }}_2\left(\mathrm{\psi }_{23}{}_{ss}-{\displaystyle \frac{{\mathrm{\theta }}_{23}}{{\mathrm{\omega }}_2{\mathrm{\omega }}_3}}\right)\hfill \\ +{\displaystyle \frac{1}{{\mathrm{\omega }}_2}}\sin {\mathrm{\omega }}_2\left(\mathrm{\psi }_{12}{}_{ss}-{\displaystyle \frac{{\mathrm{\theta }}_{12}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}}\right)+{\displaystyle \frac{1}{{\mathrm{\omega }}_3}}\sin {\mathrm{\omega }}_3\left(\mathrm{\psi }_{31}{}_{ss}-{\displaystyle \frac{{\mathrm{\theta }}_{31}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_3}}\right)\hfill \\ -{\displaystyle \frac{1}{{\mathrm{\omega }}_3}}\sin {\mathrm{\omega }}_3\left(\mathrm{\psi }_{12}{}_{ss}-{\displaystyle \frac{{\mathrm{\theta }}_{12}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}}\right)=0\hfill \end{array}$

$\begin{array}{c}\stackrel{.}{\mathrm{\psi }_{31}{}_{ss}}=-{\displaystyle \frac{1}{{\mathrm{\omega }}_3}}\sin {\mathrm{\omega }}_3\left(\mathrm{\psi }_{31}{}_{ss}-{\displaystyle \frac{{\mathrm{\theta }}_{31}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_3}}\right)\hfill \\ +{\displaystyle \frac{1}{{\mathrm{\omega }}_3}}\sin {\mathrm{\omega }}_3\left(\mathrm{\psi }_{23}{}_{ss}-{\displaystyle \frac{{\mathrm{\theta }}_{23}}{{\mathrm{\omega }}_2{\mathrm{\omega }}_3}}\right)+{\displaystyle \frac{1}{{\mathrm{\omega }}_1}}\sin {\mathrm{\omega }}_1\left(\mathrm{\psi }_{12}{}_{ss}-{\displaystyle \frac{{\mathrm{\theta }}_{12}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}}\right)\hfill \\ -{\displaystyle \frac{1}{{\mathrm{\omega }}_1}}\sin {\mathrm{\omega }}_1\left(\mathrm{\psi }_{31}{}_{ss}-{\displaystyle \frac{{\mathrm{\theta }}_{31}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_3}}\right)=0\hfill \end{array}$

The trivial and desirable solution set for the above equation is:

${\mathrm{\psi }}_{12}^{*}$=$\frac{{\mathrm{\theta }}_{12}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_2}$, ${\mathrm{\psi }}_{23}^{*}$=$\frac{{\mathrm{\theta }}_{23}}{{\mathrm{\omega }}_2{\mathrm{\omega }}_3}$,${\mathrm{\psi }}_{31}^{*}$=$\frac{{\mathrm{\theta }}_{31}}{{\mathrm{\omega }}_1{\mathrm{\omega }}_3}$ or ${\mathrm{\sigma }}_{12}^{*}$=$\mathrm{\sigma }_{23}{}^{*}$=${\mathrm{\sigma }}_{31}^{*}$=0

However, like before in eq. 13a, there are several spurious solutions for which the entire right-hand sides of eqs. 14a, 14b, 14c are zero, without each of the individual terms being zero separately. So, σ_ij can achieve any of the possible solutions satisfying the above equation. It appears that finding the spurious solutions analytically are not possible.

To verify this numerically eq. 11 was simulated for the given set of network parameters: ω₁ = 4, ω₂ = 7, ω₃ = 10, A_ij = A = 0.2, $\text{θ}=\left[\begin{array}{ccc}0& 1.3871& 1.576\\ -1.3871& 0& 0.1889\\ -1.576& -0.1889& 0\end{array}\right]$ and the initial values of ${\mathrm{\varnothing }}_i^{\prime }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 $\mathrm{\sigma }_{ij}{}_{ss}$ do not depend on θ_ij but on the natural frequencies of the oscillators (ω_i). The following Figure 7 summarizes the dependency of $\mathrm{\sigma }_{ij}{}_{ss}$ on ∅_i(0) and ∅_j(0) where it can be observed that $\mathrm{\sigma }_{12}{}_{ss}$, $\mathrm{\sigma }_{23}{}_{ss}$, and $\mathrm{\sigma }_{31}{}_{ss}$ 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 ∅₁, ∅₂ and ∅₃, $\mathrm{\sigma }_{ij}{}_{ss}$’s attained one of the solutions and for a confined subspace of ∅₁(0), ∅₂(0) and ∅₃(0) under space 0 ≤ ∅_i(0) < 2π, $\mathrm{\sigma }_{ij}{}_{ss}$’s attained the zero or desired solution $\left({\mathrm{\sigma }}_{ij}^{*}\right)$.

FIGURE 7

Figure 7. Equation 11 is simulated for 100 s for the network parameter values mentioned previously to check for what combinations of initial values of ${\mathrm{\varnothing }}_i^{\prime }s$, $\mathrm{\sigma }_{ij}{}_{ss}$’s reaches, which of the solutions of eq. 14. (A) and (B) depicts that for various initial values of ∅₁and ∅₃ in between 0 to  2π, $\mathrm{\sigma }_{31}{}_{ss}$ attains 4 possible solutions: 0.3235, −0.0036, −0.4576, −1.3652, similarly from (C-F) it can be seen that for various initial values of ∅₁, ∅₂ and ∅₂, ∅₃ in between 0 to  2π, $\mathrm{\sigma }_{12}{}_{ss}$ and $\mathrm{\sigma }_{23}{}_{ss}$ attains the 0.99,  1.465, 1.6307,  1.9426, −0.0046 and 0.0633,  0.6855, −0.6751 solutions respectively, satisfying eq. 14.

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.

${\mathrm{\tau }}_W\stackrel{.}{W_{ij}}=-W_{ij}+z_i{\left(z_j{}^{*}\right)}^{\frac{{\mathrm{\omega }}_i}{{\mathrm{\omega }}_j}}(\mathrm{15a})$

The polar coordinate representation:

${\mathrm{\tau }}_W\stackrel{.}{A_{ij}}=-A_{ij}+r_{i}r_j{}^{\frac{{\mathrm{\omega }}_i}{{\mathrm{\omega }}_j}}\cos \left({\mathrm{\omega }}_i\left(\frac{{\mathrm{\varnothing }}_i}{{\mathrm{\omega }}_i}-\frac{{\mathrm{\varnothing }}_j}{{\mathrm{\omega }}_j}-\frac{{\mathrm{\theta }}_{ij}}{{\mathrm{\omega }}_i{\mathrm{\omega }}_j}\right)\right)(\mathrm{15b})$

${\mathrm{\tau }}_W\stackrel{.}{{\mathrm{\theta }}_{ij}}=\frac{{\mathrm{\omega }}_j{r}_{i}r_j{}^{\frac{{\mathrm{\omega }}_i}{{\mathrm{\omega }}_j}}}{A_{ij}}\sin \left({\mathrm{\omega }}_i\left(\frac{{\mathrm{\varnothing }}_i}{{\mathrm{\omega }}_i}-\frac{{\mathrm{\varnothing }}_j}{{\mathrm{\omega }}_j}-\frac{{\mathrm{\theta }}_{ij}}{{\mathrm{\omega }}_i{\mathrm{\omega }}_j}\right)\right)(\mathrm{15c})$

We have limited the scope of our study only to θ_ij dynamics by assuming $\stackrel{.}{A_{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 $\stackrel{.}{A_{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) ≈ ω_it + ∅_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}^{\prime }s)$ trained according to eq. 15, θ_ij learns ${\mathrm{\varnothing }}_i\left(0\right){\mathrm{\omega }}_j-{\mathrm{\varnothing }}_j\left(0\right){\mathrm{\omega }}_i(\approx {\mathrm{\varnothing }}_i{\mathrm{\omega }}_j-{\mathrm{\varnothing }}_j\mathrm{\omega }{}_i)$.

In such a network when individual oscillators $\left(z_i^{\prime }sor{\mathrm{\varnothing }}_i^{\prime }s\right)$ are driven by complex sinusoidal forcing $(I_{ext}{}_i=I_{0i}{e}^{i\left({\mathrm{\omega }}_{i}t+{\mathrm{\phi }}_i\right)}forthe{i}^{th}oscillator)$ and the power coupling weights are trained as stated, θ_ij learns φ_iω_j−φ_jω_i as with A_ij≪1, $\stackrel{.}{A_{ij}}=0$ and I_0i≫A_ij. ${\mathrm{\varnothing }}_i^{\prime }$s are approximately same as ω_it + φ_i. The dynamics of such a network of supercritical Hopf oscillators is defined below:

$\stackrel{.}{z_i}=z_i\left(\mathrm{\mu }+i{\mathrm{\omega }}_i-{\left|z_i\right|}^2\right)+\sum _{j=1\ni j\ne i}^N{A}_{ij}{e}^{i\frac{{\mathrm{\theta }}_{ij}}{{\mathrm{\omega }}_j}}z_j{}^{\frac{{\mathrm{\omega }}_i}{{\mathrm{\omega }}_j}}+I_{ext}{}_i(16)$

$I_{ext}{}_i=I_{0i}{e}^{i\left({\mathrm{\omega }}_{i}t+{\mathrm{\phi }}_i\right)}$

Polar coordinate representation:

$\begin{array}{c}\stackrel{.}{r_i}=\left(\mathrm{\mu }-r_i{}^2\right)r_i+\sum _{j=1\ni j\ne i}^N{A}_{ij}r_j{}^{\frac{{\mathrm{\omega }}_{\text{i}}}{{\mathrm{\omega }}_{\text{j}}}}\cos {\mathrm{\omega }}_{\text{i}}\left(\frac{{\mathrm{\varnothing }}_j}{{\mathrm{\omega }}_j}-\frac{{\mathrm{\varnothing }}_i}{{\mathrm{\omega }}_i}+\frac{{\mathrm{\theta }}_{ij}}{{\text{ùù}}_{\text{i}}{\mathrm{\omega }}_{\text{j}}}\right)\hfill \\ +I_{0i}\cos \left({\mathrm{\omega }}_{i}t+{\mathrm{\phi }}_i-{\mathrm{\varnothing }}_i\right)\hfill \end{array}(\mathrm{17a})$

$\begin{array}{c}\stackrel{.}{{\mathrm{\varnothing }}_i}={\mathrm{\omega }}_i+\sum _{j=1\ni j\ne i}^N{A}_{ij}\frac{r_j{}^{\frac{{\mathrm{\omega }}_i}{{\mathrm{\omega }}_j}}}{r_i}sin{\mathrm{\omega }}_i\left(\frac{{\mathrm{\varnothing }}_j}{{\mathrm{\omega }}_j}-\frac{{\mathrm{\varnothing }}_i}{{\mathrm{\omega }}_i}+\frac{{\mathrm{\theta }}_{ij}}{{\mathrm{\omega }}_i{\mathrm{\omega }}_j}\right)\hfill \\ +\frac{I_{0i}}{r_i}\sin \left({\mathrm{\omega }}_{i}t+{\mathrm{\phi }}_i-{\mathrm{\varnothing }}_i\right)\hfill \end{array}(\mathrm{17b})$

The dynamics of such a network of Kuramoto oscillators:

$\begin{array}{c}\stackrel{.}{{\mathrm{\varnothing }}_i}={\mathrm{\omega }}_i+\sum _{j=1\ni j\ne i}^N{A}_{ij}sin{\mathrm{\omega }}_i\left(\frac{{\mathrm{\varnothing }}_j}{{\mathrm{\omega }}_j}-\frac{{\mathrm{\varnothing }}_i}{{\mathrm{\omega }}_i}+\frac{{\mathrm{\theta }}_{ij}}{{\mathrm{\omega }}_i{\mathrm{\omega }}_j}\right)\hfill \\ +I_{0i}sin\left({\mathrm{\omega }}_{i}t+{\mathrm{\phi }}_i-{\mathrm{\varnothing }}_i\right)\hfill \end{array}(18)$

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_{ext}{}_i=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).

FIGURE 8

Figure 8. When two Hopf oscillators (z₁, z₂) coupled through power coupling according to eq. 16 (for N=2 and $I_{ext}{}_i=0$), the schematic of the network is elaborated in panel (A), with the parameters ω₁ = 5, ω₂ = 10, A₁₂ = A₂₁ = 0.0001, θ₁₂(0) = 1.657, θ₂₁(0) = −1.657 and the power coupling weights are trained according to eq. 15, keeping the magnitude of the power coupling weight (A_ij) constant, ω_j times the angle of power coupling, θ_ij learns ∅_i(0)ω_j−∅_j(0)ω_i or ∅_iω_j−∅_jω_i. Here ∅₁(0) = 1.2046, ∅₂(0) = 2.7008, and τ_W = 10³. It is clear from the plot (B) that θ₁₂ reaches ∅₁ω₂−∅₂ω₁ and θ₂₁ reaches ∅₂ω₁−∅₁ω₂ w.r.t time. When two such Hopf oscillators as described in the previous simulation are forced by an external complex sinusoidal perturbation as described in eq. 16 ($I_{ext}{}_i\ne 0$) with the following parameters ω₁ = 5, ω₂ = 10, A₁₂ = A₂₁ = 0.0001, θ₁₂(0) = −2.513, θ₂₁(0) = 2.513, τ_W = 10³ and $I_{ext}{}_1=0.5e^{i\left(5t+\frac{\mathrm{\pi }}{4}\right)}$, $I_{ext}{}_2=0.5e^{i\left(10t+\frac{\mathrm{\pi }}{6}\right)}$, θ₁₂, and θ₂₁ learn φ₁ω₂−φ₂ω₁ = −5.7131 and φ₂ω₁−φ₁ω₂ = 5.7131 respectively (C). (D) It can be verified that the difference between θ₁₂ and ∅₁ω₂−∅₂ω₁ as well as θ₂₁ and ∅₂ω₁−∅₁ω₂ becomes zero w.r.t time.

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_{ext}{}_i$ is the complex sinusoidal input signal driving i^th oscillator, then the normalized phase difference among them is:

$\mathrm{\psi }_{ij}{}^I={\displaystyle \frac{angle\left(I_{ext}{}_i\right)}{{\mathrm{\omega }}_{0i}}}-{\displaystyle \frac{angle\left(I_{ext}{}_j\right)}{{\mathrm{\omega }}_{0j}}}$

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).

$\stackrel{.}{z_i}=z_i\left(\mathrm{\mu }+i{\mathrm{\omega }}_i-{\left|z_i\right|}^2\right)+\sum _{j=1\ni j\ne i}^N{A}_{ij}{e}^{i\frac{{\mathrm{\theta }}_{ij}}{{\mathrm{\omega }}_j}}z_j{}^{\frac{{\mathrm{\omega }}_i}{{\mathrm{\omega }}_j}}+\mathrm{\epsilon }I_{ext}{}_i(\mathrm{19a})$

$I_{ext}{}_i=I_{0i}{e}^{i\left({\mathrm{\omega }}_{0i}t+{\mathrm{\phi }}_i\right)}$

${\mathrm{\tau }}_W\stackrel{.}{W_{ij}}=-W_{ij}+z_i{\left(z_j{}^{*}\right)}^{\frac{{\mathrm{\omega }}_i}{{\mathrm{\omega }}_j}}(\mathrm{19b})$

$\stackrel{.}{{\mathrm{\omega }}_i}=-\mathrm{\epsilon }{I}_{0i}\sin \left({\mathrm{\varnothing }}_i-{\mathrm{\omega }}_{0i}t-{\mathrm{\phi }}_i\right)(\mathrm{19c})$

FIGURE 9

Figure 9. Equations 19a-c are simulated for the following set of parameters:ω₀₁ = 20, ω₀₂ = 30, ω₁(0) = 30, ω₂(0) = 40, A₁₂ = A₂₁ = 0.0001, θ₁₂(0) = −1.7884, θ₂₁(0) = 1.7884, τ_W = 10³, ε = 0.9 and $I_{ext}{}_1=e^{i\left(20t+\frac{\mathrm{\pi }}{4}\right)}$, $I_{ext}{}_2=e^{i\left(30t+\frac{\mathrm{\pi }}{6}\right)}$. It can be verified that θ_ij learns φ_jω_0i−φ_iω_0j after the ${\mathrm{\omega }}_i^{\prime }$s learns the corresponding ${\mathrm{\omega }}_{0i}^{\prime }$s.

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.

To this end, we construct a reservoir of Hopf oscillators (Figure 10A). It consists of a network of Hopf oscillators with trainable lateral connections (eq. 15) and trainable natural frequencies (eq. 20b). In addition, there exists a linear weight stage between the oscillators and the output layer. The natural frequencies of the oscillators adapt themselves to the nearest significant component in the input signal. The lateral connections involving power coupling, encode the (normalized) phase relationships among the oscillators. The output weights represent the amplitudes of the oscillatory components corresponding to the oscillators. A similar network architecture and dynamics for adaptable central pattern generator (aCPG) has been proposed by Righetti et al. (2005), the detailed comparison of which is discussed in Section “Discussion.”

FIGURE 10

Figure 10. The schematic of the network where (A) a reservoir of N number of adaptive Hopf oscillators similar to the network proposed by 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_pi is the predicted output signal at the i^th output node. Complex feed-forward weights $W_{ij}^{f\prime }s$ connecting j^th oscillator to the i^th output node, are trained in batch mode according to eq. 23.

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.

$\stackrel{.}{z_i}=z_i\left(\mathrm{\mu }+i{\mathrm{\omega }}_i-{\mathrm{\beta }}_1{\left|z_i\right|}^2\right)+\sum _{j=1\ni j\ne i}^N{W}_{ij}z_j{}^{\frac{{\mathrm{\omega }}_i}{{\mathrm{\omega }}_j}}+\mathrm{\epsilon }e\left(t\right)(\mathrm{20a})$

$W_{ij}=A_{ij}{e}^{i\frac{{\mathrm{\theta }}_{ij}}{{\mathrm{\omega }}_j}}$

$\stackrel{.}{{\mathrm{\omega }}_i}=-{\mathrm{\eta }}_{\mathrm{\omega }}e\left(t\right)\sin {\mathrm{\varnothing }}_i(\mathrm{20b})$

$\stackrel{.}{{\mathrm{\alpha }}_i}={\mathrm{\eta }}_{\mathrm{\alpha }}e\left(t\right)r_i\cos {\mathrm{\varnothing }}_i(\mathrm{20c})$

${\mathrm{\tau }}_W\stackrel{.}{W_{ij}}=-W_{ij}+z_i{\left(z_j{}^{*}\right)}^{\frac{{\mathrm{\omega }}_i}{{\mathrm{\omega }}_j}}\mathit{\hspace{1em}}\ni \stackrel{.}{A_{ij}}=0(\mathrm{20d})$

$e\left(t\right)=D\left(t\right)-P(t)(20e)$

$P\left(t\right)={\displaystyle \sum _{i=1}^N}{\mathrm{\alpha }}_i\cos {\mathrm{\varnothing }}_i$

Let D(t) be the teaching time-series signal with a finite number of frequencies ($D\left(t\right)={\sum }_{i=1}^N{a}_i\cos \left({\mathrm{\omega }}_{i}t+{\mathrm{\phi }}_i\right)$ where φ_i be the phase offset associated with i^th frequency component), z_i = r_ie^i∅_i, ${\mathrm{\alpha }}_i^{\prime }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.

FIGURE 11

Figure 11. Equations 20a-e is simulated with the following parameters $D(t)=2\cos \left(4t+\frac{\mathrm{\pi }}{2}\right)+1.5\cos \left(8t+\frac{\mathrm{\pi }}{5}\right)+1.8\cos \left(12t+\frac{\mathrm{\pi }}{12}\right)$, A_ij = 10⁻⁵, η_ω = 0.1, η_α = 10⁻⁴, τ_W = 10⁴, ε = 0.5, N = 3, dt = 10⁻³ sec. (A) 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.

We propose a similar network for Fourier decomposition of complex time series signal with a finite number of frequencies (assuming $D\left(t\right)={\sum }_{i=1}^N{a}_i{e}^{i\left({\mathrm{\omega }}_{i}t+{\mathrm{\phi }}_i\right)}$), 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.

$\stackrel{.}{z_i}=z_i\left(\mathrm{\mu }+i{\mathrm{\omega }}_i-{\mathrm{\beta }}_1{\left|z_i\right|}^2\right)+\sum _{j=1\ni j\ne i}^N{W}_{ij}z_j{}^{\frac{{\mathrm{\omega }}_i}{{\mathrm{\omega }}_j}}+\mathrm{\epsilon }e\left(t\right)(\mathrm{21a})$

$W_{ij}=A_{ij}{e}^{i\frac{{\mathrm{\theta }}_{ij}}{{\mathrm{\omega }}_j}}$

$\stackrel{.}{{\mathrm{\omega }}_i}=-{\mathrm{\eta }}_{\mathrm{\omega }}\left(real\left(e\left(t\right)\right)\sin {\mathrm{\varnothing }}_i-img\left(e\left(t\right)\right)\cos {\mathrm{\varnothing }}_i\right)(\mathrm{21b})$

$\stackrel{.}{{\mathrm{\alpha }}_i}={\mathrm{\eta }}_{\mathrm{\alpha }}\left(real\left(e\left(t\right)\right)r_i\cos {\mathrm{\varnothing }}_i+img\left(e\left(t\right)\right)r_i\sin {\mathrm{\varnothing }}_i\right)(\mathrm{21c})$

${\mathrm{\tau }}_W\stackrel{.}{W_{ij}}=-W_{ij}+z_i{\left(z_j{}^{*}\right)}^{\frac{{\mathrm{\omega }}_i}{{\mathrm{\omega }}_j}}\ni \stackrel{.}{A_{ij}}=0(\mathrm{21d})$

$e\left(t\right)=D\left(t\right)-P(t)(21e)$

$P\left(t\right)={\displaystyle \sum _{i=1}^N}{\mathrm{\alpha }}_i{z}_i$

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).

FIGURE 12

Figure 12. Equations 21a-e are simulated for a teaching signal constituting three frequency components with the parameters $D\left(t\right)=2e^{i\left(4t+\frac{\mathrm{\pi }}{2}\right)}+1.5e^{i\left(8t+\frac{\mathrm{\pi }}{5}\right)}+1.8e^{i\left(12t+\frac{\mathrm{\pi }}{12}\right)},A_{ij}={10}^{-5},{\mathrm{\eta }}_{\mathrm{\omega }}=0.1,{\mathrm{\eta }}_{\mathrm{\alpha }}={10}^{-4},{\mathrm{\tau }}_W={10}^4,\mathrm{\epsilon }=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.

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 (${\mathrm{\omega }}_i^{\prime }$s), the real feed-forward weights (${\mathrm{\alpha }}_i^{\prime }s$) as well as the power coupling weights $\left(W_{ij}^{\prime }s\right)$ 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 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.

$W_{ij}^f=K_{ij}{e}^{i{\mathrm{\zeta }}_{ij}}(\mathrm{22a})$

$Y_p{}_i\left(t\right)=real\left(\sum _{j=1}^N{W}_{ij}^f{e}^{i{\mathrm{\varnothing }}_j}\right)(\mathrm{22b})$

$\mathrm{△}{K}_{ij}=\left(-1\right){\mathrm{\eta }}_K\sum _t\left(Y_d{}_i\left(t\right)-Y_p{}_i\left(t\right)\right)\cos \left({\mathrm{\varnothing }}_j(t)+{\mathrm{\zeta }}_{ij}\right)(\mathrm{22c})$

$\mathrm{△}{\mathrm{\zeta }}_{ij}=\left(-1\right){\mathrm{\eta }}_{\mathrm{\zeta }}\sum _t\left(Y_d{}_i\left(t\right)-Y_p{}_i\left(t\right)\right)\left(-K_{ij}\sin \left({\mathrm{\varnothing }}_j\left(t\right)+{\mathrm{\zeta }}_{ij}\right)\right)(\mathrm{22d})$

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 $\frac{{\mathrm{\omega }}_i}{{\mathrm{\omega }}_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.

FIGURE 13

Figure 13. (A) In the first phase of training the oscillatory reservoir network (Figure 10A) with three oscillators is trained using the following teaching signal and learning parameters: D(t) = A₁cos⁡(ω₁t + ∅₁) + A₂cos⁡(ω₂t + ∅₂) + A₃cos⁡(ω₂t + ∅₂), A_ij = 10⁻¹, η_w = 10⁻¹, η_α = 10⁻⁴, τ_w = 10⁴, ϵ = 0.5, N=3, dt = 10⁻³sec, n_epoch = 15. The network produced and the original teaching signal and their magnitude spectrum respectively after the last epoch (15) of training. (B) After the second phase of training, the network (Figure 10B) is able to reconstruct the M=5 output signals (Y_di 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⁻⁵, η_ζ = 10⁻⁶, no. of epochs for the batch mode of learning = 10000.

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.

FIGURE 14

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⁻¹ for |j−i| < 5 and A_ij = 0 for |j−i|≥5η_ω = 0.1, η_α = 10⁻⁴, τ_W = 10⁴, ε = 0.5, N = 100, dt = 10⁻³ 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⁻⁵, η_ζ = 10⁻⁶, no. of epochs for the batch mode of learning = 1000.

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 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 long-range 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 real-valued counterparts to the complex domain. Complex feed-forward 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 multi-oscillator networks with real coupling.

Network models of complex-valued oscillators have been described before. For example, subsequent studies by Hoppensteadt and 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 (ω₁ ≠ ω₂) 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ω₁ = nω₂, 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∅₁−n∅₂ for Ω₁₂ = mω₁−nω₂ = 0, and increases/decreases linearly or rotates for Ω₁₂ ≠ 0 at steady state. However, neither of the study extended the result to the case of interacting oscillators with very different natural frequencies (|mω₁−nω₂| > ϵ, 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 (∅_i−∅_j, m∅_i−n∅_j, $\frac{{\mathrm{\varnothing }}_i}{{\mathrm{\omega }}_i}-\frac{{\mathrm{\varnothing }}_j}{{\mathrm{\omega }}_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).

TABLE 1

Table 1. The brief comparison among the various coupling strategies between a pair of complex Hopf oscillators.

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, $\frac{{\mathrm{\phi }}_i}{{\mathrm{\omega }}_i}-\frac{{\mathrm{\phi }}_j}{\mathrm{\omega }{}_j}$.

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 ${\mathrm{\alpha }}_i{}^{\prime }\text{s}$ and ${\mathrm{\omega }}_i{}^{\prime }\text{s}$, when the network is reset (i.e., $r_i{}^{\prime }\text{s}$ and ${\mathrm{\varnothing }}_i{}^{\prime }\text{s}$ are reset) there is no external forcing to drive the phase offset of each oscillator to the desired phase offsets $\left({\mathrm{\phi }}_i^{\prime }\text{s}\right)$. 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 $\left(\frac{{\mathrm{\varnothing }}_i}{{\mathrm{\omega }}_i}-\frac{{\mathrm{\varnothing }}_0}{{\mathrm{\omega }}_0}\right)$. 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ω₀, ω₀ 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 above-mentioned 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−ω_it) 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 $\mathrm{\sigma }_{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.

TABLE 2

Table 2. This table compares various features related to network architecture, scope of trainability, retrievability, type of teaching signals the model can learn; of the two adaptive central pattern generator (aCPG) models proposed by Righetti et al. (2005) with our proposed aCPG model.

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⁡(ω₀t + ∅_x(t))). The three characteristic features, envelop (Z_x(t)), constant mean frequency (ω₀) and phase (∅_x(t)) respectively are dependent on network parameters. ω₀ 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{\mathrm{\theta }}_{\mathbf{p}}}$

where, θ_p is a N-dimensional vector, N being the no of oscillators. The power coupling coefficients are assigned as:

$W_{ij}=a{e}^{i\frac{{\mathrm{\Theta }}_{ij}}{{\mathrm{\omega }}_j}}$

Where,

$e^{i\mathrm{\Theta }}=\sum _p{x}_{p}x_p{}^{*}$

$x_p{}^{*}$ is the complex conjugate of the transpose of x_p. However, the network retains the phase relationship among the oscillators defined by the combined phase, Θ, instead of the individual phase relationships defined by various patterns. The model equations in the current form do not seem to support storage of multiple patterns within the oscillator layer. However, with slightly modified equations in the complex domain, it is possible to achieve multiple pattern storage (Chakravarthy and Ghosh, 1996; Chakravarthy, 2009). The proposed network architectures seem to be robust to the noise inherent to the system dynamics as well as external noise, corrupting the input signal. A brief study on the performance of key network architectures w.r.t the noise power is given in the Supplementary Figures 1-5.

Future Perspective

In the future, we plan to extend the proposed model to create a whole class of deep oscillatory networks. The network will have an input stage that serves as an encoder that performs a Fourier-like decomposition of the input time-series signals. The network model of Section “A Network for Reconstructing a Signal by a Fourier-Like Decomposition” will play the role of this encoder. The hidden layers will consist of oscillators operating at a range of frequencies. The output layer will be a decoder that converts the oscillatory outputs of the last hidden layer into the output time series. Another interesting proposed study is to develop the network model of Section “A Network for Reconstructing a Signal by a Fourier-Like Decomposition” into a model of the tonotopic map. Electrophysiological recordings from the bat’s auditory cortex revealed a map of frequencies (a tonotopic map). Kohonen had proposed a model of the tonotopic map using a self-organizing map (Kohonen, 1998). However, this model represents frequency as an explicitly available parameter, without actually modeling the responses of the neurons to oscillatory inputs (tones). We propose that by organizing the oscillators in the reservoir model of Section “A Network for Reconstructing a Signal by a Fourier-Like Decomposition” in a 2D array with neighborhood connections, it is possible to produce a biologically more feasible tonotopic map model.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

The simulations were done mostly by DB with the support of SP. DB wrote the main text. VC contributed to providing the key ideas, conducting initial simulations to test the hypothesis, editing the manuscript drafts, and providing insight into structure. All authors contributed to the article and approved the submitted version.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We acknowledge the support from undergraduate student from the Mechanical Engineering Department of IIT Madras, Asit Tarsode for his insight on the dependency of steady state solutions on the initial conditions. We also acknowledge the constant and much needed support from senior lab mate Vignayanandam R Muddapu, and other lab members. This manuscript has been released as a pre-print at bioRxiv (Biswas et al., 2020) for feedback from fellow neuroscientists.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2021.551111/full#supplementary-material

References

Auer, P., Burgsteiner, H., and Maass, W. (2008). A learning rule for very simple universal approximators consisting of a single layer of perceptrons $. Neural Netw. 21, 786–795. doi: 10.1016/j.neunet.2007.12.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Bi, G. Q., and Poo, M. M. (1998). Synaptic modifications in cultured hippocampal neurons: dependence on spike timing, synaptic strength, and postsynaptic cell type. J. Neurosci. 18, 10464–10472. doi: 10.1523/jneurosci.18-24-10464.1998

PubMed Abstract | CrossRef Full Text | Google Scholar

Biswas, D., Sooryakiran, P., and Chakravarthy, V. S. (2020). A complex-valued oscillatory neural network for storage and retrieval of multichannel electroencephalogram signals. bioRxiv [Preprint]. doi: 10.1101/2020.03.28.013292

CrossRef Full Text | Google Scholar

Breakspear, M., Heitmann, S., and Daffertshofer, A. (2010). Generative models of cortical oscillations: neurobiological implications of the Kuramoto model. Front. Hum. Neurosci. 4:190. doi: 10.3389/fnhum.2010.00190

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunel, N., and Wang, X. J. (2003). What determines the frequency of fast network oscillations with irregular neural discharges? I. Synaptic dynamics and excitation-inhibition balance. J. Neurophysiol. 90, 415–430. doi: 10.1152/jn.01095.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsáki, G., Anastassiou, C. A., and Koch, C. (2012). The origin of extracellular fields and currents — EEG, ECoG, LFP and spikes. Nat. Rev. Neurosci. 13, 407–420. doi: 10.1038/nrn3241

PubMed Abstract | CrossRef Full Text | Google Scholar

Campbell, S., and Wang, D. (1994). Synchronization and desynchronization in a network of locally coupled synchronization and desynchronization in a network of locally coupled Wilson-Cowan oscillators. IEEE Trans. Neural Netw. 7, 541–554. doi: 10.1109/72.501714

CrossRef Full Text | Google Scholar

Chakravarthy, S. V., and Ghosh, J. (1996). A complex-valued associative memory for storing patterns as oscillatory states. Biol. Cybern. 75, 229–238. doi: 10.1007/s004220050290

CrossRef Full Text | Google Scholar

Chakravarthy, V. S. (2009). “A complex-valued hopfield neural network: dynamics and applications,” in Complex-Valued Neural Networks: Utilizing High-Dimensional Parameters, ed. T. Nitta, (Hershey, PA: Information Science Reference), 79–103. doi: 10.4018/978-1-60566-214-5.ch004

CrossRef Full Text | Google Scholar

Chakravarthy, V. S. (2020). Encoding temporal relations using complex – valued synapses. ACCS 1, 1–3.

Google Scholar

Craik, A., and He, Y. (2019). Deep learning for electroencephalogram (EEG) classification tasks?: a review. J. Neural Eng. 16, 031001. doi: 10.1088/1741-2552/ab0ab5

PubMed Abstract | CrossRef Full Text | Google Scholar

Damoiseaux, J. S., and Greicius, M. D. (2009). Greater than the sum of its parts: a review of studies combining structural connectivity and resting-state functional connectivity. Brain Struct. Funct. 213, 525–533. doi: 10.1007/s00429-009-0208-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Destexhe, A., Mainen, Z. F., and Sejnowski, T. J. (1994). Synthesis of models for excitable membranes, synaptic transmission and neuromodulation using a common kinetic formalism. J. Comput. Neurosci. 1, 195–230. doi: 10.1007/BF00961734

PubMed Abstract | CrossRef Full Text | Google Scholar

Draguhn, A. (2004). Neuronal oscillations in cortical networks. Science 304, 1926–1930. doi: 10.1126/science.1099745

PubMed Abstract | CrossRef Full Text | Google Scholar

Ermentrout, B., and Kopell, N. (1991). Multiple pulse interactions and averaging in systems of coupled neural oscillators^∗. J. Math. Biol. 29, 195–217. doi: 10.1007/BF00160535

CrossRef Full Text | Google Scholar

Ermentrout, G. B., and Kopell, N. (1989). Mathematical Modelling of Central Pattern Generators. Cambridge, MA: Academic Press Limited, doi: 10.1016/b978-0-12-287960-9.50013-7

CrossRef Full Text | Google Scholar

Frasconi, P., Gori, M., and Soda, G. (1995). Recurrent neural networks and prior knowledge for sequence processing?: a constrained nondeterministic approach. Knowledge Based Syst. 8, 313–332. doi: 10.1016/0950-7051(96)81916-2

CrossRef Full Text | Google Scholar

Fujii, H., Ito, H., Aihara, K., Ichinose, N., and Tsukada, M. (1996). Dynamical cell assembly hypothesis – theoretical possibility of spatio-temporal coding in the cortex. Neural Netw. 9, 1303–1350. doi: 10.1016/S0893-6080(96)00054-8

CrossRef Full Text | Google Scholar

Ghosh-dastidar, S., and Lichtenstein, A. G. (2009). Review article spiking neural networks. Int. J. 19, 295–308.

Google Scholar

Grandchamp, R., Braboszcz, C., and Delorme, A. (2014). Oculometric variations during mind wandering. Front. Psychol. 5:31. doi: 10.3389/fpsyg.2014.00031

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoff, A., Santos, J. V., Manchein, C., and Albuquerque, H. A. (2014). Numerical bifurcation analysis of two coupled FitzHugh-Nagumo oscillators. Eur. Phys. J. B. 87:151. doi: 10.1140/epjb/e2014-50170-9

CrossRef Full Text | Google Scholar

Hopfield, J. J. (1982). Neural networks and physical systems with emergent collective computational abilities. Proc. Natl. Acad. Sci. U.S.A. 79, 2554–2558. doi: 10.1073/pnas.79.8.2554

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoppensteadt, F. C., and Izhikevich, E. M. (1996). Biological Cybernetics Synaptic organizations and dynamical properties of weakly connected neural oscillators. II. Learning phase information. Biol. Cybern. 135, 129–135. doi: 10.1007/s004220050280

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoppensteadt, F. C., and Izhikevich, E. M. (1997). Weakly Connected Neural Networks. New York, NY: Springer-Verlag.

Google Scholar

Hoppensteadt, F. C., and Izhikevich, E. M. (2000). Pattern recognition via synchronization in phase-locked loop neural networks. IEEE Trans. Neural Netw. 11, 734–738. doi: 10.1109/72.846744

CrossRef Full Text | Google Scholar

Hyafil, A., Giraud, A. L., Fontolan, L., and Gutkin, B. (2015). Neural cross-frequency coupling: connecting architectures, mechanisms, and functions. Trends Neurosci. 38, 725–740. doi: 10.1016/j.tins.2015.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Izhikevich, E. M. (2003). Simple model of spiking neurons. IEEE Trans. Neural Netw. 14, 1569–1572. doi: 10.1109/TNN.2003.820440

PubMed Abstract | CrossRef Full Text | Google Scholar

Izhikevich, E. M. (2004). Which model to use for cortical spiking neurons? IEEE Trans. Neural Netw. 15, 1063–1070. doi: 10.1109/TNN.2004.832719

PubMed Abstract | CrossRef Full Text | Google Scholar

Izhikevich, E. M. (2007). Dynamical Systems in Neuroscience?: The Geometry of Excitability and Bursting. Cambridge, MA: The MIT Press.

Google Scholar

Jansen, B. H., and Rit, V. G. (1995). Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biol. Cybern. 73, 357–366. doi: 10.1007/BF00199471

PubMed Abstract | CrossRef Full Text | Google Scholar

Joglekar, M. R., Mejias, J. F., Yang, G. R., and Wang, X. J. (2018). Inter-areal balanced amplification enhances signal propagation in a large-scale circuit model of the primate cortex. Neuron 98, 222–234.e8. doi: 10.1016/j.neuron.2018.02.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, J. C., and Large, E. W. (2015). Signal processing in periodically forced gradient frequency neural networks. Front. Comput. Neurosci. 9:152. doi: 10.3389/fncom.2015.00152

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, J. C., and Large, E. W. (2021). Multifrequency Hebbian plasticity in coupled neural oscillators. Biol. Cybern. 115, 43–57. doi: 10.1007/s00422-020-00854-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Kohonen, T. (1998). The self-organizing map. Neurocomputing 21, 1–6.

Google Scholar

Lawrence, S., Giles, C. L., Tsoi, A. C., and Back, A. D. (1997). Face recognition: a convolutional neural-network approach. IEEE Trans. Neural Netw. 8, 98–113. doi: 10.1109/72.554195

CrossRef Full Text | Google Scholar

Lippmann, R. P. (1990). Review of neural networks for speech recognition. Read. Speech Recognit. 38, 374–392. doi: 10.1016/b978-0-08-051584-7.50036-x

CrossRef Full Text | Google Scholar

Low, L. A., Reinhall, P. G., Storti, D. W., and Goldman, E. B. (2005). Coupled van der Pol oscillators as a simplified model for generation of neural patterns for jellyfish locomotion. Struct. Control Heal. Monit. 13, 417–429. doi: 10.1002/stc.133

CrossRef Full Text | Google Scholar

Maass, W. (1997a). Fast sigmoidal networks via spiking neurons. Neural Comput. 9, 279–304. doi: 10.1162/neco.1997.9.2.279

PubMed Abstract | CrossRef Full Text | Google Scholar

Maass, W. (1997b). Networks of spiking neurons: the third generation of neural network models. Neural Netw. 10, 1659–1671. doi: 10.1016/S0893-6080(97)00011-7

CrossRef Full Text | Google Scholar

Mejias, J. F., Murray, J. D., Kennedy, H., and Wang, X. J. (2016). Feedforward and feedback frequency-dependent interactions in a large-scale laminar network of the primate cortex. Sci. Adv. 2:e1601335. doi: 10.1126/sciadv.1601335

PubMed Abstract | CrossRef Full Text | Google Scholar

Pikovsky, A., Rosenblum, M., and Kurths, J. (2001). Synchronization A Universal Concept in Nonlinear Sciences. Cambridge: Cambridge University Press.

Google Scholar

Powanwe, A. S., and Longtin, A. (2019). Determinants of brain rhythm burst statistics. Sci. Rep. 9:18335. doi: 10.1038/s41598-019-54444-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Righetti, L., Buchli, J., and Ijspeert, A. J. (2005). “From dynamic hebbian learning for oscillators to adaptive central pattern generators,” in Proceedings of 3rd International Symposium on Adaptive Motion in Animals and Machines – AMAM 2005, (Ilmenau: Verlag ISLE), 1–7.

Google Scholar

Righetti, L., Buchli, J., and Ijspeert, A. J. (2006). Dynamic Hebbian learning in adaptive frequency oscillators. Phys. D Nonlinear Phenom. 216, 269–281. doi: 10.1016/j.physd.2006.02.009

CrossRef Full Text | Google Scholar

Roy, Y., Banville, H., and Albuquerque, I. (2019). Deep learning-based electroencephalography analysis?: a systematic review. J. Neural Eng. 16:051001. doi: 10.1088/1741-2552/ab260c

PubMed Abstract | CrossRef Full Text | Google Scholar

Ruck, D. W., Rogers, S. K., and Kabrisky, M. (1990). Feature selection using a multilayer perceptron. J. Neural Netw. Comput. 2, 40–48.

Google Scholar

Sadilek, M., and Thurner, S. (2015). Physiologically motivated multiplex Kuramoto model describes phase diagram of cortical activity. Sci. Rep. 5:10015. doi: 10.1038/srep10015

PubMed Abstract | CrossRef Full Text | Google Scholar

Saha, A., and Feudel, U. (2017). Intermittency in delay-coupled FitzHugh – Nagumo oscillators and loss of phase synchrony as its precursor. Indian Acad. Sci. Conf. Ser. 1:1. doi: 10.29195/iascs.01.01.0010

CrossRef Full Text | Google Scholar

Schmidhuber, J. (2015). Deep learning in neural networks?: an overview. Neural Netw. 61, 85–117. doi: 10.1016/j.neunet.2014.09.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, Z. W., Shi, Z. Z., Liu, X., and Shi, Z. P. (2008). A computational model for feature binding. Sci. China C Life Sci. 51, 470–478. doi: 10.1007/s11427-008-0063-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Stefanescu, R. A., and Jirsa, V. K. (2008). A low dimensional description of globally coupled heterogeneous neural networks of excitatory and inhibitory neurons. PLoS Comput. Biol. 4:e1000219. doi: 10.1371/journal.pcbi.1000219

PubMed Abstract | CrossRef Full Text | Google Scholar

Toral, R., Mirasso, C., and Gunton, J. (2003). System size coherence resonance in coupled FitzHugh-Nagumo models. Europhys. Lett. 61:162. doi: 10.1209/epl/i2003-00207-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Trappenberg, T. P. (2003). Continuous Attractor Neural Networks Biographical Sketch?: Continuous Attractor Neural Networks Abstract. (Hershey, PA: IGI Global, 701 E. Chocolate Ave), 1–30.

Google Scholar

Wang, D., and Terman, D. (1997). Image segmentation based on oscillatory correlation. Neural Comput. 9, 805–836. doi: 10.1162/neco.1997.9.4.805

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, D. L., and Terman, D. (1995). Locally excitatory globally inhibitory oscillator networks. IEEE Trans. Neural Netw. 6, 283–286. doi: 10.1109/72.363423

CrossRef Full Text | Google Scholar

Wilson, H. R., and Cowan, J. D. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J. 12, 1–24. doi: 10.1016/S0006-3495(72)86068-5

CrossRef Full Text | Google Scholar

Wirkus, S., and Rand, R. (2002). The dynamics of two coupled van der pol oscillators. Nonlinear Dyn. 30, 205–221.

Google Scholar

Wong, K. F., and Wang, X. J. (2006). A recurrent network mechanism of time integration in perceptual decisions. J. Neurosci. 26, 1314–1328. doi: 10.1523/JNEUROSCI.3733-05.2006

PubMed Abstract | CrossRef Full Text | Google Scholar