Figure 1: Phase of oscillation (denoted by \(\vartheta\) in the rest of the article) of the FitzHugh-Nagumo model with I=0.5. The phase here is measured in units of time. The zero-phase point \(x_0\) is chosen to correspond to the peak of the potential (the peak of spike).

Coupled oscillators interact via mutual adjustment of their amplitudes and phases. When coupling is weak, amplitudes are relatively constant and the interactions could be described by phase models.

Phase of oscillation

Many physical, chemical, and biological systems can produce rhythmic oscillations (Winfree 2001), which can be represented mathematically by a nonlinear dynamical system
\[
x' = f(x)
\]
having a periodic orbit \(\gamma\ .\) Let \(x_0\) be an arbitrary point on \(\gamma\ ,\) then any other point on the periodic orbit can be characterized by the time, \(\vartheta\ ,\) since the last passing of \(x_0\ ;\) see Figure 1. The variable \(\vartheta\) is called phase of oscillation, and it is bounded by the period of oscillation \(T\ .\) The phase is often normalized by \(T\) or \(T/2\pi\ ,\) so that it is bounded by \(1\) or \(2\pi\ ,\) respectively.

The phase of oscillation can also be defined outside \(\gamma\) using the notion of isochrons. The change of variables \(x(t) = \gamma(\vartheta(t))\) transforms the nonlinear system in a neighborhood of \(\gamma\) into an equivalent but simpler phase model
\[
\vartheta' = 1.
\]
Such a change of variables removes the amplitude but saves the phase of oscillation (Ermentrout 1986). It is often convenient to assume that the phase \(\vartheta\) is defined on the unit circle \(S^1\ .\)

The phase of oscillation could also be defined for chaotic systems (Pikovsky et al. 2001) using the observation that many chaotic attractors, as in Rossler oscillator, look like smeared limit cycles.

Weak forcing

The same change of variables transforms a weakly forced oscillator
\[
x' = f(x) + \varepsilon s(t)
\]
into the phase model of the form
\[
\vartheta' = 1 + \varepsilon Q(\vartheta) \cdot s(t),
\]
where, the term \(\varepsilon s(t)\) denotes a weak time-dependent (and possibly \(x\)-dependent) input, e.g., from other oscillators in a network; the dot, "\(\cdot\)", denotes the scalar (dot) product of two vectors; The function \(Q(\vartheta)\ ,\) illustrated in the Figure 2, is called linear response function, sensitivity function, or infinitesimal PRC. This impulse response function, kernel, or Green's function can then be convolved with the actual input received by each oscillator in order to compute the total phase resetting of the oscillator received over one cycle of the network oscillation. \(Q(\vartheta)\) satisfies three equivalent conditions (Izhikevich 2007):

Winfree\[Q(\vartheta)\] is a normalized phase response curve (PRC) to infinitesimal pulsed perturbations. That is, one measures PRC of the oscillator \(x'=f(x)\) by perturbing each component of state vector \(x\) with brief pulses of small amplitude with area of the pulse \(A\ ,\) and then takes \(Q=\)PRC\(/A\) in the limit \(A\rightarrow 0 \ .\)

Kuramoto\[Q(\vartheta) = \]grad\(\Theta(x)\ ,\) where \(\Theta(x)\) is the isochron function defined in a neighborhood of the periodic orbit \(\gamma\ .\) That is, one starts from every point \(x\) in a neighborhood of \(\gamma\) and determines its asymptotic phase, \(\Theta(x)\ ,\) relative to the phase of the solution starting with \(x_0\ .\)

Malkin\[Q(\vartheta)\] is the solution to the adjoint problem \(
dQ/d\vartheta=-\{Df(\gamma(\vartheta))\}^\top Q\;, \) with the normalization condition \(Q(\vartheta) \cdot f(\gamma(\vartheta))=1\) for any \(\vartheta\ .\) That is, one determines the Jacobian matrix \(Df\) along the periodic orbit and then solves, usually numerically, the adjoint problem.

Malkin's condition, though least intuitive, is the most useful in applications.

Examples of reduction

The infinitesimal PRC function \(Q(\vartheta)\) can be found analytically in a few simple cases.

The implicit assumption of weak coupling is that the relative position (phase) of the oscillators changes slowly with respect to their motion around the limit cycle (absolute phase). This implies slow convergence to a steady state phase locking. If the coupling is not sufficiently weak but is pulsatile in nature, the methods for pulse-coupled oscillators can be utilized, otherwise there are no general methods.

Phase model

Introducing phase deviation variables \(\vartheta_i = t + \varphi_i\ ,\) one can transform the system above into the form
\[
\varphi_i' = \varepsilon \sum_{j=1}^n h_{ij}(t + \varphi_i, \ t + \varphi_j).
\]
This system can be averaged to the phase model
\[
\varphi_i' = \varepsilon \sum_{j=1}^n H_{ij}(\varphi_j-\varphi_i)\;,
\]
where each function
\[
H_{ij}(\chi) = \lim_{T\rightarrow \infty} \frac{1}{T} \int_0^T h_{ij}(t, \ t +\chi)\, dt
\]
describes the interaction between oscillators. This function is just a constant unless the oscillators have nearly resonant periods, i.e., the ratio \(T_i/T_j\) is \(\varepsilon\)-close to a low-order rational number \(p/q\) (\(p+q\) is small). Since the dynamics of two coupled non-resonant oscillators is described by an uncoupled phase model (\(H=\)const), such oscillators cannot phase lock. That is, the phase of one of them cannot change the phase of the other one even on the long time scale of order \(1/\varepsilon\ .\)

Computational neuroscience provides an important application of phase models. In this case, the state variables \(x_i\) and \(x_j\) describe activities of the post-synaptic (forced) and pre-synaptic (forcing) periodically firing neurons, and the function \(g_{ij}\) describes the time course of synaptic input. The phase variables \(\varphi_i\) and \(\varphi_j\) describe the timings of firings of the neurons, and the function \(H_{ij}\) describes normalized phase resetting curve (Netoff et al 2005).

Figure 3: Examples of the connection function H; modified from Izhikevich (2007). The units on the y-axis are "phase/time".

where
\[
\omega = \omega_2 - \omega_1
\]
and
\(
H(\chi) = H_{21}(-\chi)-H_{12}(\chi),
\)
is the frequency mismatch and the anti-symmetric part of the coupling, respectively, illustrated in the Figure 3, dashed curves. A stableequilibrium of this system corresponds to a stable limit cycle of the phase model.

All equilibria of this system are solutions to \(H(\chi) = -\omega\ ,\) provided that \(\omega\) is small enough. Geometrically, the equilibria are intersections of the horizontal line \(-\omega\) with the graph of \(H\ .\) They are stable if the slope of the graph is negative at the intersection, i.e., \(H'(\chi)<0\ .\) If oscillators are identical, then \(H(\chi)\) is an odd function (i.e., \(H(-\chi)=-H(\chi)\)), and \(\chi=0\) and \(\chi=\pi\) are always equilibria, possibly unstable, corresponding to the in-phase and anti-phase synchronized solutions. The in-phase synchronization of coupled oscillators in the figure is stable because the slope of \(H\) (dashed curves) is negative at \(\chi=0\ .\) The max and min values of the function \(H\) determine the tolerance of the network to the frequency mismatch \(\omega\ ,\) since there are no equilibria outside this range.

Chains of oscillators

The behavior of chains of phase models is considerably more complex than that of pairs, even for nearest neighbor coupling. The reason for this is that when coupling is local, oscillators at the ends get different inputs from those in the middle so that phase locking may not even exist. However, in a large class of models, chains can be analyzed either by direct calculation or by letting the size of the chains tend to infinity. In the former case, Cohen et al. (1982) examined a linear chain of nearest neighbor oscillators with a frequency gradient:
\[
\theta_i' = \omega_i + \sin (\theta_{i+1}-\theta_i) + \sin(\theta_{i-1}-\theta_i).
\]
As long as the differences in the frequencies are small enough, there will be a phase-locked solution. Interestingly, if the length of the chain is \( N \) and the frequency gradient is linear with slope \(b\) then \(b = O(1/N)\) as \(N\to\infty\ .\) That is, nearest neighbor chains can support very small gradients when the coupling is sinusoidal (and, in fact, any odd periodic function). However, if the coupling function contains any even components (that is, replace \( sin\theta \) with \( \sin (\theta+\beta) \ ,\) then frequency gradients as that are \( O(1) \) can be supported in nearest neighbor chains of coupled phase oscillators. Kopell & Ermentrout (1986,1990) derived a set of continuum equations from which general phase-locked solutions could be found.

Networks of oscillators that are not arranged in a ring have "boundaries" that can lead to patterns of phase difference that look like waves. If the coupling is isotropic, the waves take the form of one-dimensional target waves, either originating at the center and propagating symmetrically to the edges, or starting at the edges and propagating to the center. For example, an array of nearest neighbor-coupled oscillators of the form
\[
\theta_i' = \sin(\theta_{i+1}-\theta_i - \alpha) + \sin(\theta_{i-1}-\theta_i -\alpha) \quad i=0,\ldots,N
\]
with no "i-1" term for oscillator 0 and no "i+1" term for oscillator N, will generate a wave with phase which are roughly of the form \(
\theta_i = \alpha |N/2-i|. \)
With anisotropic coupling, the waves are almost straight lines, \( \theta_i = \pm \alpha i. \)
These results were proven in a series of papers by Ermentrout and Kopell.

Linear arrays of oscillators

Now consider a network of \(n>2\) weakly all to all coupled oscillators. To determine the existence and stability of synchronized states in the network, we need to study equilibria of the corresponding phase model
\[
\varphi_i' = \varepsilon \omega_i + \varepsilon \sum_{j\neq i}^n H_{ij}(\varphi_j-\varphi_i).
\]
Existence of one equilibrium of the phase model above implies the existence of the entire circular family of equilibria, since translation of all \(\varphi_i\) by a constant phase shift does not change the phase differences \(\varphi_i-\varphi_j\) and hence the form of the phase model. This family corresponds to a periodic orbit, on which all oscillators have equal frequencies and constant phase shifts, i.e., they are synchronized, possibly out-of-phase.

Vector \(\phi=(\phi_1,\ldots,\phi_n)\) is an equilibrium when
\[
0 = \omega_i + \sum_{j\neq1}^n H_{ij}(\phi_j-\phi_i)
\] for all \(i\ .\)
It is stable when all eigenvalues of the linearization matrix (Jacobian) at \(\phi\) have negative real parts, except one zero eigenvalue corresponding to the eigenvector along the circular family of equilibria (\(\phi\) plus a phase shift is a solution too since the phase shifts \(\phi_j-\phi_i\) are not affected).

In general, determining the stability of equilibria is a difficult problem. Ermentrout (1992) found a simple sufficient condition. If

\(a_{ij} = H_{ij}'(\phi_j-\phi_i) \geq 0\ ,\) and

the directed graph defined by the matrix \(a = (a_{ij})\) is connected, (i.e., each oscillator is influenced, possibly indirectly, by every other oscillator),

then the equilibrium \(\phi\) is neutrally stable, and the corresponding limit cycle \(x(t+\phi)\) of the phase model is asymptotically stable.

Another sufficient condition was found by Hoppensteadt and Izhikevich (1997). If the phase model satisfies

\(\omega_1=\cdots = \omega_n = \omega\) (identical frequencies)

\(H_{ij}(-\chi) = - H_{ji}(\chi)\) (pair-wise odd coupling)

for all \(i\) and \(j\ ,\) then the network dynamics converge to a limit cycle. On the cycle, all oscillators have equal frequencies \(1+\varepsilon\omega\) and constant phase deviations. The proof follows from the observation that the phase model is a gradient system in a rotating coordinate system.

2D Arrays of oscillators

Two-dimensional arrays provide a much richer class of dynamics. Consider the two-dimensional analogues of the one-dimensional chain of nearest-neighbor coupled oscillators:
\[
\theta_{i,j}' = H_N(\theta_{i+1,j}-\theta_{i,j}) + H_S(\theta_{i-1,j}-\theta_{i,j})
+ H_E(\theta_{i,j+1}-\theta_{i,j}) + H_W(\theta_{i,j-1}-\theta_{i,j})
\]
It can be shown that the patterns of phase-locked activity are of the form
\[
\theta_{i,j} = \Psi_i + \Phi_j
\]
where \( \Psi_i,\Phi_j \) are solutions to the one-dimensional chains using \( H_{N,S} \) or \( H_{E,W} \) respectively.