# A Phase Model with Large Time Delayed Coupling111This research is supported
in part by the Natural Sciences and Engineering Research Council of Canada.

Isam Al-Darabsah222Department of Applied Mathematics, University of Waterloo,
Waterloo, ON, N2L 3G1, Canada. ,333Email<EMAIL_ADDRESS>Sue Ann
Campbell$\star$,444Email<EMAIL_ADDRESS>
###### Abstract

We consider two identical oscillators with weak, time delayed coupling. We
start with a general system of delay differential equations then reduce it to
a phase model. With the assumption of large time delay, the resulting phase
model has an explicit delay and phase shift in the argument of the phases and
connection function, respectively. Using the phase model, we prove that for
any type of oscillators and any coupling, the in-phase and anti-phase phase-
locked solutions always exist and give conditions for their stability. We show
that for small delay these solutions are unique, but with large enough delay
multiple solutions of each type with different frequencies may occur. We give
conditions for the existence and stability of other types of phase-locked
solutions. We discuss the various bifurcations that can occur in the phase
model as the time delay is varied. The results of the phase model analysis are
applied to Morris-Lecar oscillators with diffusive coupling and compared with
numerical studies of the full system of delay differential equations. We also
consider the case of small time delay and compare the results with the
existing ones in the literature.

Keywords: Coupled oscillators $\cdot$ Large time delay $\cdot$ Synchronization
$\cdot$ Phase-locking

##  1 Introduction

Coupled oscillator models have been used to study different aspects of
biology, chemistry and engineering, for example chemical waves [1], flashing
of fireflies [2], laser arrays [3, 4], power system networks [5], neural
networks [6, 7, 8, 9], movement of a slime mold [10], and coupled predator-
prey systems [11, 12]. Time delays in the connections between the oscillators
are inescapable due to the time for a signal to propagate from one element to
the other. Many of these systems exhibit phase-locking behaviour, i.e., all
the oscillators have similar waveforms and frequencies, but with some fixed
phase difference between different oscillators. To study the existence and
stability of such phase-locked solutions and how they are related to the time
delay and other parameters, one must formulate a model for the system. We
discuss two approaches below.

One approach to study connected networks of oscillators is through phase
models [13]. In these models, each oscillator is represented only by its phase
along its limit cycle, with amplitude variation neglected [14, 15, 16]. Phase
models take the general form [17, 14]:

$\frac{d\theta_{i}}{d\xi}=\Omega_{i}+H_{i}\left(\theta_{1}(\xi),\ldots,\theta_{n}(\xi)\right),\quad
i=1,\ldots,n,$ (1)

where $\theta_{i}\in[0,2\pi)$ is the phase of the $i^{\rm th}$ oscillator,
$\Omega_{i}>0$ the natural frequency and $H_{i}$ are the connection functions.
Motivated by the famous Kuramoto model [1], in the literature the functions
$H_{i}$ often take the form:

$H_{i}\left(\theta_{1}(\xi),\ldots,\theta_{n}(\xi)\right)=\sum_{j=1}^{n}K_{ij}H\left(\theta_{j}(\xi)-\theta_{i}(\xi)\right),\quad
i=1,\ldots,n,$ (2)

where $K_{ij}$ is the adjacency matrix of an unweighted network [15, 18]. In
the original Kuramoto model [1] the function $H$ in (2) is the sine function.
Usually, transmission time delay is introduced as an explicit delay in the
argument of the phases [18, 19, 20, 21, 22, 23]:

$\frac{d\theta_{i}}{d\xi}=\Omega_{i}+\sum_{j=1}^{n}K_{ij}H\left(\theta_{j}(\xi-\tau)-\theta_{i}(\xi)\right),\quad
i=1,\ldots,n.$ (3)

Most studies of this model focus only on synchronization [18] or use
simplifications such as $H(\cdot)=\sin(\cdot)$ [24, 19, 20, 21, 23, 22] or
$n=2$ [19, 23, 22].

Other authors introduce additional processes into system (3). For instance, in
[19], the dynamic behavior of coupled oscillators with time delayed
interaction under a pinning force is studied. In [21, 22], the authors study
time delayed phase models with $H=\sin(\cdot)$ and random noise forcing.
Finally, a phase shift is sometimes included in the model of a network of
connected oscillators to represent the temporal distance between the
oscillators. In general, the phase shift between two oscillators $\alpha_{ij}$
is incorporated in the phase model as, see e.g., [25, 24, 26],

$\frac{d\theta_{i}}{d\xi}=\Omega_{i}+\sum_{j=1}^{n}K_{ij}H\left(\theta_{j}(\xi)-\theta_{i}(\xi)-\alpha_{ij}\right),\quad
i=1,\ldots,n.$ (4)

In the case where $H(\cdot)=\sin(\cdot)$ this model is called the Kuramoto-
Sakaguchi model[26]. In fact, there is a relation between such phase shifts
and the transmission time delay. In [27, 28], the authors have shown how the
model with delay and the model with the phase shift are linked. We will review
the details of this link later in this section.

Models of coupled oscillators are also formulated as physically or biological
derived differential equations [11, 12, 17]. These models are of the form

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{d\mathbf{X}_{i}}{d\rho}={\mathbf{F}}_{i}\left(\mathbf{X}_{i}(\rho)\right)+\epsilon{\mathbf{G}}_{i}\big{(}\mathbf{X}_{1}(\rho),\ldots,\mathbf{X}_{i}(\rho),\ldots,\mathbf{X}_{n}(\rho)\big{)},\quad
i=1,\ldots,n,\quad\mathbf{X}_{i}\in\mathbb{R}^{m},}$ (5)

and are such that when $\epsilon=0$ the dynamical system of each uncoupled
oscillator has an exponentially asymptotically stable $T_{i}-$periodic limit
cycle with corresponding (natural) frequency $\Omega_{i}$. In these models,
$\mathbf{X}_{i}$ represents the state of the $i^{th}$ oscillator of the
system, ${\mathbf{G}}_{i}$ are the coupling functions and $\epsilon>0$ is the
coupling strength [29, 14, 27, 30]. Note that $\mathbf{X}_{i}$ is a vector of
dimension at least $2$, but can be high dimensional. For example, in a
pendulum model $\mathbf{X}_{i}$ represents the position and velocity of the
$i^{th}$ pendulum, while in a neural model $\mathbf{X}_{i}$ represents the
voltage and gating variables of the $i^{th}$ neuron.

If the coupling is weak, $0<\epsilon\ll 1$, then the theory of weakly coupled
oscillators can be used to connect the physical model (5) to a phase model [8,
29, 30, 31, 32] . More precisely, the dynamics of each oscillator in the
network can be rigorously reduced to a single equation that indicates how the
phase of the oscillator changes in time [14, 27, 16]. One form of weakly
coupled oscillator theory is Malkin’s Theorem where the connection functions
in the phase model are determined explicitly in terms of ${\mathbf{G}}_{i}$
and the limit cycles of the uncoupled system, (5) with $\epsilon=0$. Let
$\varphi_{i}(t)\in\mathbb{S}^{1}$ be the phase deviation of the $i^{\rm th}$
oscillator of (5), i.e., the change in the phase due to the coupling. It then
follows from Malkin’s Theorem (see e.g., [14, Theorem 9.2]) that the dynamics
of (5) can be described by the phase deviation model:

$\displaystyle\frac{d\varphi_{i}}{d{t}}$
$\displaystyle=H_{i}\big{(}\varphi_{1}\left(t\right)-\varphi_{i}\left(t\right),\ldots,\varphi_{n}\left(t\right)-\varphi_{i}\left(t\right)\big{)}+\mathcal{O}(\epsilon),\quad
i=1,\ldots,n,$ (6)

where $H_{i}$ are the phase interaction functions and the variable
$t:=\epsilon\rho$ represents slow time because the phase deviations
$\varphi_{i}$ are slow variables. The references [30, 14, 16] provide other
forms of the theory and give further references. We also refer the reader to
the recent articles [33, 34, 35] for an overview of various numerical and
analytical techniques for phase reduction. In [27], Izhikevich generalizes
Malkin’s theorem to weakly connected oscillators with fixed delay, $\tau$, in
their interaction:

$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{d\mathbf{X}_{i}}{d\rho}={\mathbf{F}}_{i}\left(\mathbf{X}_{i}(\rho)\right)+\epsilon{\mathbf{G}}_{i}\big{(}\mathbf{X}_{1}(\rho-\tau),\ldots,\mathbf{X}_{i}(\rho-\tau),\ldots\mathbf{X}_{n}(\rho-\tau)\big{)},\quad
i=1,\ldots,n,\quad\mathbf{X}_{i}\in\mathbb{R}^{m}}$ (7)

where all uncoupled oscillators have nearly identical natural frequencies.
Assuming the natural frequency is $1$, Izhikevich shows that the phase
deviation model corresponding to (7) is

$\displaystyle\frac{d\varphi_{i}}{d{t}}$
$\displaystyle=H_{i}\big{(}\varphi_{1}\left(t-\eta\right)-\varphi_{i}\left(t\right)-\zeta,\ldots,\varphi_{n}\left(t-\eta\right)-\varphi_{i}\left(t\right)-\zeta\big{)}+\mathcal{O}(\epsilon),\quad
i=1,\ldots,n,$ (8)

where $\eta:=\epsilon\tau$ and $\zeta:=\tau\mod 2\pi$. The functions $H_{i}$
are still defined explicitly in terms of ${\mathbf{G}}_{i}$ and the uncoupled
limit cycle in (7). It is clear that the time delay $\tau$ enters the phase
model (8) as both an explicit delay, $\eta$, and a phase shift, $\zeta$. The
major result that Izhikevich proved in [27] is that if the delay $\tau$ in (7)
satisfies $\epsilon\tau=\mathcal{O}(1)$ (large delay), then the explicit delay
occurs in the phase model (8). However, when the delay satisfies
$\tau=\mathcal{O}(1)$ with respect to $\epsilon$ (small delay), no delay
appears in the argument of the phases. Hence, (8) becomes:

$\displaystyle\frac{d\varphi_{i}}{d{t}}$
$\displaystyle=H_{i}\big{(}\varphi_{1}\left(t\right)-\varphi_{i}\left(t\right)-\zeta,\ldots,\varphi_{n}\left(t\right)-\varphi_{i}\left(t\right)-\zeta\big{)}+\mathcal{O}(\epsilon),\quad
i=1,\ldots,n.$ (9)

We refer the reader to the review article [24] and the references therein for
different scenarios where large or small delay appears in-phase models.

In this article we focus on physical models with the following particular form

$\frac{{d{{\mathbf{X}}_{i}}}}{{d\rho}}={\mathbf{F}}({{\mathbf{X}}_{i}}(\rho))+\epsilon\sum\limits_{j=1}^{n}{{K_{ij}}}{\mathbf{G}}({{\mathbf{X}}_{i}}(\rho),{{\mathbf{X}}_{j}}(\rho-\tau)),\quad
i=1,\ldots n,\ {{\mathbf{X}}_{i}}\in{\mathbb{R}^{m}}$ (10)

where $K_{ii}=0$. This represent the following modelling assumptions. The
oscillators are identical. The coupling occurs pairwise between the
oscillators and there is no coupling from an oscillator to itself. The
coupling to the $i^{th}$ oscillator occurs close to that oscillator, so the
time delay represents the time it takes for information to travel from the
$j^{th}$ oscillator to the $i^{th}$ oscillator. Models with such structure
occur in models of biological systems [8, 11].

Assuming the uncoupled oscillators in (10) have a natural frequency $\Omega$
and the $K_{ij}=\mathcal{O}(1)$ with respect to $\epsilon$, we show in the
appendix that the approach of [27] can be applied to yield

$\frac{d\varphi_{i}}{dt}=\frac{1}{\Omega}\sum_{j=1}^{n}K_{ij}H(\varphi_{j}(t-\eta)-\varphi_{i}(t)-\zeta)+\mathcal{O}(\epsilon)$
(11)

where $\eta:=\epsilon\Omega\tau$ and $\zeta:=\Omega\tau\mod 2\pi$, in the case
of large delay, i.e., when $\epsilon\Omega\tau=\mathcal{O}(1)$. In the case of
small delay (11) becomes

$\frac{{d{\varphi_{i}}(t)}}{{dt}}=\frac{1}{{\Omega}}\sum\limits_{j=1}^{n}{{K_{ij}}}H\left({{\varphi_{j}}(t)-{\varphi_{i}}(t)-\Omega\tau}\right)+\mathcal{O}(\epsilon).$
(12)

To see how the phase deviation model relates to the standard phase model, note
that the phase of oscillations $\theta_{i}$ in (10) have the form:

$\theta_{i}(\xi)=\Omega\xi+\varphi_{i}(t),\quad i=1,\ldots,n,$ (13)

where $t=\epsilon\Omega\xi$, see [27, 14]. Notice that the natural frequency
of each uncoupled oscillator in (13) is $\Omega$. Then,

$\displaystyle\frac{d\theta_{i}}{d\xi}=\Omega+\epsilon\Omega\frac{d\varphi_{i}}{dt}=\Omega+\epsilon\sum_{j=1}^{n}K_{ij}H\left(\theta_{j}(\xi-\tau)-\theta_{i}(\xi)\right)+\mathcal{O}(\epsilon^{2}).$
(14)

Similarly, when the time delay is small, we have

$\displaystyle\frac{d\theta_{i}}{d\xi}=\Omega+\epsilon\sum_{j=1}^{n}K_{ij}H\left(\theta_{j}(\xi)-\theta_{i}(\xi)-\zeta\right)+\mathcal{O}(\epsilon^{2}).$
(15)

Thus in the phase model formulation, the coupling strength parameter
$\epsilon$ explicitly appears in front of the connection function $h$.
Regarding the dynamics, it follows from (13) that

$\theta_{i+1}-\theta_{i}=\varphi_{i+1}-\varphi_{i},\quad i=1,\ldots,{n-1}$

i.e., phase-locked solutions are the same as phase deviation locked solutions
[14]. The existence and stability of phase-locked solutions of system (15) has
been studied in the case of two oscillators [36, 24] and many oscillators with
structured coupling [17, 24, 37].

The goals in this paper are twofold. First, the majority of studies of coupled
oscillators with large delays have been done in the context of isolated phase
models, often with just sine function coupling. Thus we will revisit and
extend this analysis in the case where the phase model is explicitly connected
to a physical differential equation model and the function $H$ is general. In
particular, we will show that the multiple stable phase-locked solutions of
the same type may occur even when the coupling is weak. Second, note that the
small delay phase deviation model (12) is a system of ordinary differential
equations, while the large delay model (10) is a delay differential equation
model. Thus the spectrum of Floquet multipliers of a periodic solution is
finite for the former and countably infinite for the latter. Nevertheless,
several studies have verified numerically that the model (12) gives an
accurate description of existence and stability of phase-locked periodic
solutions of (5) in the case of weak coupling and small delay [36, 17]. Here
we will show why this is the case. In particular we will show how the
solutions of system (10) reduce to those of system (12) if the delay is small.
In this article, we will focus on (10) when $n=2$ as this is enough to
illustrate our main points.

The paper is organized as follows. In the next section, we reduce the model of
two weakly connected oscillators with large time delay to a phase model, and
study the existence of phase-locked solutions. In Section 3, we give a
complete description of the stability criteria for all phase-locked solutions
and describe the potential bifurcations that can occur in the system. Then we
compare our results with the stability criteria in [36] when the time delay is
small. In Section 4, we consider a particular application to Morris-Lecar
oscillators with diffusive coupling. Numerically, we derive the corresponding
phase model, calculate the phase-locked solutions, determine their stability
and explore the existence of bifurcations. We also compare prediction of the
phase model and solutions of the full model. Finally, we examine the behaviour
when the time delay is small. In Section 5, we discuss our results.

##  2 Phase Model

Consider the system of ODEs

$\frac{d\mathbf{X}_{i}}{d\rho}={\mathbf{F}}({\mathbf{X}_{i}}(\rho))\quad
i=1,2,\quad\mathbf{X}_{i}\in\mathbb{R}^{n}.$ (16)

Assume that the system (16) admits an exponentially asymptotically stable
periodic orbit given by $\mathbf{X}=\tilde{\mathbf{X}}(\rho)$ with natural
frequency $\Omega$, $0\leq\rho\leq T=2\pi/\Omega$.

Next, consider a weakly connected system of two identical coupled oscillators
of the form (16) with time delayed coupling:

$\displaystyle\frac{d\mathbf{X}_{1}}{d{\rho}}$
$\displaystyle=\mathbf{F}\left(\mathbf{X}_{1}(\rho)\right)+\epsilon\mathbf{G}\left(\mathbf{X}_{1}(\rho),\mathbf{X}_{2}(\rho-\tau);\epsilon\right),$
(17) $\displaystyle\frac{d\mathbf{X}_{2}}{d\rho}$
$\displaystyle=\mathbf{F}\left(\mathbf{X}_{2}(\rho)\right)+\epsilon\mathbf{G}\left(\mathbf{X}_{2}(\rho),\mathbf{X}_{1}(\rho-\tau);\epsilon\right),$

where $\mathbf{G}:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$
describes the coupling between the two oscillators and $\epsilon$ is the
coupling strength. Assume that $\epsilon$ is sufficiently small and
$\eta:=\epsilon\Omega\tau=\mathcal{O}(1)$. Let ${t}=\epsilon\rho$ be slow time
and $\varphi_{i}(t)\in\mathbb{S}^{1}$ be the phase deviation from the natural
oscillation $\hat{X}(\rho)$, $\rho\geq 0$. Then, by applying weakly coupled
oscillator theory with delayed interactions in [27],
$(\varphi_{1},\varphi_{2})^{T}\in\mathbb{T}^{2}$ is a solution to

$\displaystyle\frac{d\varphi_{1}}{d{t}}$
$\displaystyle=\frac{1}{\Omega}H(\varphi_{2}(t-\eta)-\varphi_{1}(t)-\Omega\tau)+\mathcal{O}(\epsilon),$
(18) $\displaystyle\frac{d\varphi_{2}}{d{t}}$
$\displaystyle=\frac{1}{\Omega}H(\varphi_{1}(t-\eta)-\varphi_{2}(t)-\Omega\tau)+\mathcal{O}(\epsilon),$

where $H$ is a $2\pi-$periodic function defined by

$H(\phi)=\frac{1}{{2\pi}}\int\limits_{0}^{2\pi}{\hat{\mathbf{Z}}{{(\rho)}^{T}}\mathbf{G}\left({\hat{\mathbf{X}}(\rho),\hat{\mathbf{X}}(\rho+\phi)}\right)}d\rho.$
(19)

Here $\hat{\mathbf{Z}}{(\rho)}$ is the unique nontrivial $2\pi-$periodic
solution to the adjoint linear system

$\frac{{d\hat{\mathbf{Z}}}}{{d\rho}}=-{\left[{D\mathbf{F}\left({\hat{\mathbf{X}}(\rho)}\right)}\right]^{T}}\hat{\mathbf{Z}}$

satisfying the normalization condition

$\frac{1}{{2\pi}}\int\limits_{0}^{2\pi}{\hat{\mathbf{Z}}(\rho)\cdot}\mathbf{F}\left({\hat{\mathbf{X}}(\rho)}\right)d\rho=1.$

The derivation of system (18) from (17) follows from the Appendix with $n=2$
and $K_{12}=K_{21}=1$.

Dropping the terms $\mathcal{O}(\epsilon)$ in (18), we obtain the phase
deviation model:

$\displaystyle\frac{d\varphi_{1}}{d{t}}$
$\displaystyle=\frac{1}{\Omega}H(\varphi_{2}(t-\eta)-\varphi_{1}(t)-\Omega\tau),$
(20) $\displaystyle\frac{d\varphi_{2}}{d{t}}$
$\displaystyle=\frac{1}{\Omega}H(\varphi_{1}(t-\eta)-\varphi_{2}(t)-\Omega\tau).$

For simplicity, in the rest of the paper we will refer to (20) as the phase
model instead of the phase deviation model.

We study the dynamics of the model (20) by exploring phase locking in (20),
that is, solutions of (20) such that $\varphi_{2}-\varphi_{1}=\text{constant}$
[14]. We suppose that

$\varphi_{1}(t)=\omega t\qquad\text{and}\qquad\varphi_{2}(t)=\omega t+\psi$
(21)

where $\omega$ is the frequency deviation of the oscillator and $\psi$ is the
natural phase difference [14]. Substituting (21) into (20) leads to

$\displaystyle\omega-\frac{1}{\Omega}H(\psi-\omega\eta-\Omega\tau)$
$\displaystyle=0,$ (22)
$\displaystyle\omega-\frac{1}{\Omega}H(-\psi-\omega\eta-\Omega\tau)$
$\displaystyle=0.$

We rewrite this as

$\displaystyle F(\omega,\psi)=0=F(\omega,-\psi)$ (23)

where

$\displaystyle
F(\omega,\cdot):=\omega-\frac{1}{\Omega}H(\cdot-\omega\eta-\Omega\tau).$ (24)

In this article, we are interested in exploring how the solutions ($\psi$ and
$\omega$) of (22) vary with $\tau$ when the coupling strength ($\epsilon$) and
frequency ($\Omega$) are fixed. Note that, we need only to investigate $\psi$
in $[0,2\pi)$, due to the $2\pi$ periodicity of $H$, and
$\omega\in\mathbb{R}$.

First, by subtracting the equations of (22), we obtain

$H(\psi-\omega\eta-\Omega\tau)-H(-\psi-\omega\eta-\Omega\tau)=0.$ (25)

Since $H$ is $2\pi-$periodic function, equation (25) always has the solutions
$\psi=0,\pi$. The corresponding frequency deviation is determined from the
equation

$F(\omega,0)=\omega-\frac{1}{\Omega}H(-\omega\eta-\Omega\tau)=0$ (26)

when $\psi=0$ and

$F(\omega,\pi)=\omega-\frac{1}{\Omega}H(\pi-\omega\eta-\Omega\tau)=0$ (27)

when $\psi=\pi$.

Equations (26) and (27) are guaranteed to have a least one solution due to the
continuity and $2\pi$ periodicity of $H$. In fact, if $\tau$ is sufficiently
large, they may have multiple solutions. To see this, recall that
$\eta=\epsilon\omega\tau$ and note that

$\displaystyle F_{\omega}(\omega,0)=1+\epsilon\tau
H^{\prime}(-\omega\epsilon\Omega\tau-\Omega\tau),$ (28)

where $F_{\omega}$ is the partial derivative of $F$ with respect to $\omega$.
If there exists $\overline{\omega}$ such that $F(\overline{\omega},0)=0$ and
$F_{\omega}(\overline{\omega},0)<0$ then (26) has more than one solution.
Similar arguments apply to equation (27). This may be possible if $\tau$ is
sufficiently large.

######  Remark 2.1.

The solutions $\psi^{*}=0$ and $\psi^{*}=\pi$ of (25) correspond to in-phase
and anti-phase periodic solutions of the original model (17), respectively. By
in-phase solution we mean both oscillators reach their highest peak at the
same time, whereas an anti-phase solution means one oscillator reaches its
highest peak one half-period after the other oscillator. Examples of these
solutions are given in Figure 1.

Figure 1: Illustrations of the phase-locked dynamics of model (17).

In fact, system (17) could have other phase-locked solutions (neither in-phase
nor anti-phase) corresponding to the solutions $\psi$ of (22) such that
$\psi\notin\\{0,\pi\\}$. As in [38], we will refer to these solutions of (17)
as out-of-phase solutions. Let $(\omega^{*},\psi^{*})$ be a solution of (22)
at $\tau=\tau^{*}$ such that $\psi^{*}\notin\\{0,\pi\\}$. Then $\omega^{*}$
and $\psi^{*}$ satisfy (23), that is, $(\omega^{*},\psi^{*})$ is an
intersection point of the contours $F(\omega,\psi)=0$ and $F(\omega,-\psi)=0$
in the $\omega\psi-$plane. Suppose that $\psi^{*}\in(0,\pi)$ with a
corresponding $\omega^{*}$ are solutions to (22) at $\tau=\tau^{*}$, then

$\displaystyle\frac{1}{\Omega}H(\psi^{*}-2\pi-\omega^{*}\eta^{*}-\Omega\tau^{*})$
$\displaystyle=\frac{1}{\Omega}H(\psi^{*}-\omega^{*}\eta^{*}-\Omega\tau^{*})=\omega^{*}$

and

$\displaystyle\frac{1}{\Omega}H(2\pi-\psi^{*}-\omega^{*}\eta^{*}-\Omega\tau^{*})$
$\displaystyle=\frac{1}{\Omega}H(-\psi^{*}-\omega^{*}\eta^{*}-\Omega\tau^{*})=\omega^{*}$

due to the periodicity of $H$. Thus, $2\pi-\psi^{*}$ is also a solution in
(22) with corresponding $\omega^{*}$.

This leads to the following.

######  Proposition 2.1 (Existence of phase-locked solutions).

For any interaction function $H$ and any values of $\Omega$, $\epsilon$ and
$\tau$ the phase model (20) has the solutions $\psi^{*}=0$ and $\psi^{*}=\pi$
with corresponding frequency deviations determined by (26) and (27),
respectively. If $\psi^{*}\in(0,\pi)$ with corresponding $\omega^{*}$ are
solutions to (22) at $\tau=\tau^{*}$ then so is $2\pi-\psi^{*}$ with
$\omega^{*}$, i.e., solutions come in pairs.

##  3 Stability

In this section, we discuss the linear stability of the solutions (21) of
(20). The linearization of (20) about the solution (21) is

$\displaystyle\frac{{d{u_{1}}}}{{dt}}$
$\displaystyle=-a{u_{1}}(t)+a{u_{2}}(t-\eta),$ (29)
$\displaystyle\frac{{d{u_{2}}}}{{dt}}$
$\displaystyle=-b{u_{2}}(t)+b{u_{1}}(t-\eta),$

where

$a=\frac{1}{\Omega}H^{\prime}(\psi-\omega\eta-\Omega\tau)\qquad\text{and}\qquad
b=\frac{1}{\Omega}H^{\prime}(-\psi-\omega\eta-\Omega\tau).$ (30)

In (30), $H^{\prime}$ represents the derivative of $H$ with respect to its
argument. It is useful for our analysis to scale time so the delay becomes
one. Applying the scaling

$\eta s=t,U_{1}(s)=u_{1}(t),\ U_{2}(s)=u_{2}(t),\ $

results in

$\displaystyle\frac{{d{U_{1}}}}{{ds}}$ $\displaystyle=-\eta a{U_{1}}(s)+\eta
a{U_{2}}(s-1),$ (31) $\displaystyle\frac{{d{U_{2}}}}{{ds}}$
$\displaystyle=-\eta b{U_{2}}(s)+\eta b{U_{1}}(s-1).$

It follows that the corresponding characteristic equation is

$\Delta(\lambda;\eta)={\lambda^{2}}+\eta(a+b)\lambda+\eta^{2}ab-\eta^{2}ab{e^{-2\lambda}}=0.$
(32)

In the following we study the distribution of roots of this equation.

######  Proposition 3.1.

Assume $ab=0$. Then $\Delta(\lambda;\eta)$ has:

  1. i.

One positive root and one zero root when $a+b<0$;

  2. ii.

Two zero roots when $a+b=0$;

  3. iii.

One negative root and one zero root when $a+b>0$.

###### Proof.

The characteristic equation in this case reduces to

${\lambda^{2}}+\eta(a+b)\lambda=0.$

The result follows. ∎

######  Proposition 3.2.

$\Delta(\lambda;\eta)$ has a positive real root when one of the following
holds.

  1. i.

$ab>0$ and $a+b<0$;

  2. ii.

$ab<0$ and $a+b\leq 0$;

  3. iii.

$ab<0$, $a+b>0$ and $a+b+2\eta ab<0$.

###### Proof.

Define

$f(\lambda)=(\lambda+\eta a)(\lambda+\eta b)\quad\text{and}\quad
g(\lambda)=\eta^{2}abe^{-2\lambda}.$ (33)

Then $f(0)=g(0)=\eta^{2}ab$ and

$\Delta(\lambda;\eta)=0\quad\iff\quad f(\lambda)=g(\lambda).$ (34)

  1. i.

It follows from $ab>0$ and $a+b<0$ that $a<0$ and $b<0$. Since (32) is
symmetric in $a$ and $b$, without loss of generality, we may assume $b<a<0$.
Note that $f(-\eta b)=0<g(-\eta b)$. Since $f$ is positive and increasing for
$\lambda>-\eta b>0$ and $g$ is positive and decreasing for $\lambda>0$, there
exists $\lambda^{*}>-\eta b$ such that $f(\lambda^{*})=g(\lambda^{*})$, see
Figure 2(a).

  2. ii.

Assume $a>0$ and $b<0$. When $a+b<0$, $f$ is decreasing for
$\lambda\in\left(0,-\frac{a+b}{2}\eta\right)$ and is increasing for
$\lambda>-\frac{a+b}{2}\eta$. Further, $g$ increases for $\lambda>0$ and
$\lim_{\lambda\to\infty}g(\lambda)=0$, thus there exists
$\lambda^{*}\in\left(-\frac{a+b}{2}\eta,-\eta b\right)$ such that
$f(\lambda^{*})=g(\lambda^{*})$, see Figure 2(b). When $a+b=0$,
$f(0)=g(0)=\eta^{2}ab$, $f^{\prime}(0)=0<g^{\prime}(0)$ and $f$ is increasing
for $\lambda>0$. Thus, with the same arguments, $\lambda^{*}$ lies in
$\left(0,-\eta b\right)$.

  3. iii.

Assume $a>0$ and $b<0$. In this case $f$ and $g$ are increasing for
$\lambda>0$ and $g<0$ for $\lambda\geq 0$. Since
$f^{\prime}(0)=\eta(a+b)<-2\eta^{2}ab=g^{\prime}(0)$, then there exists
$\lambda^{*}\in(0,-\eta b)$ such that $f(\lambda^{*})=g(\lambda^{*})$, see
Figure 2(c).

∎

(a) $ab>0$ and $a+b<0$.

(b) $ab<0$ and $a+b<0$.

(c) $ab<0$ and $a+b>0$.

Figure 2: Positive real roots in $\Delta(\lambda;\eta)=0$.

######  Proposition 3.3.

When $ab>0$ and $a+b>0$, $\Delta(\lambda;\eta)$ has no roots with positive
real part.

###### Proof.

Since $ab>0$ and $a+b>0$, we have $a>0$ and $b>0$. Assume there is a root
$\lambda^{*}=x+iy$ of $\Delta(\lambda;\eta)=0$ with $x>0$. Then, it follows
from (33) and (34) that

$|f(\lambda^{*})|=|g(\lambda^{*})|.$ (35)

Notice that, due to the positivity of $x$ we get

$\left|{f\left({{\lambda^{*}}}\right)}\right|=\sqrt{{{(x+\eta
a)}^{2}}+{y^{2}}}\sqrt{{{(x+\eta b)}^{2}}+{y^{2}}}>{\eta^{2}}ab$

and

$\left|{g\left({{\lambda^{*}}}\right)}\right|={\eta^{2}}ab{e^{-2x}}<{\eta^{2}}ab.$

Hence, $|f(\lambda^{*})|>|g(\lambda^{*})|$, which contradicts (35). Thus, all
roots of $\Delta(\lambda;\eta)=0$ have nonpositive real parts when $ab>0$ and
$a+b>0$. ∎

######  Proposition 3.4.

$\lambda=0$ is a root of (32) for any $\eta$. If
$\eta\neq\eta^{*}:=-\frac{a+b}{2ab}$ then $\lambda=0$ is a simple root.
Otherwise, it is a double root. The double multiplicity of $\lambda=0$ occurs
only in the following cases.

  1. i.

$ab>0$ and $a+b<0$;

  2. ii.

$ab<0$ and $a+b>0$.

###### Proof.

It is clear that $\Delta(0;\eta)=0$ and
$\Delta^{\prime}(0;\eta)=\eta\left(a+b+2ab\eta\right)$ where ′ is the
derivative with respect to $\lambda$. If $\eta\neq\eta^{*}$ then
$\Delta^{\prime}(0;\eta)\neq 0$, and hence $\lambda=0$ is a simple root. When
$\eta=\eta^{*}$, we have $\Delta^{\prime}(0;\eta^{*})=0$ and

$\Delta^{\prime\prime}(0;\eta^{*})=-\frac{a^{2}+b^{2}}{ab}\neq 0.$

Thus, $\lambda=0$ has double multiplicity.

It is clear that $\eta^{*}$ exists if and only if

$-\frac{a+b}{2ab}>0\iff\\{ab>0\ \text{and}\ a+b<0\\}\ \text{or}\ \\{ab<0\
\text{and}\ a+b>0\\}.$

∎

######  Proposition 3.5.

When $ab<0$, $a+b>0$ and $a+b+2\eta ab\geq 0$, $\Delta(\lambda;\eta)$ has no
roots with positive real part.

###### Proof.

Note that the characteristic equation (32) can be written as

${\Delta}(\lambda;\eta)=\lambda^{2}+\eta(a+b)\lambda+\eta^{2}ab\int_{0}^{2}\lambda
e^{-u\lambda}du=0.$

Suppose that ${\Delta}(\lambda;\eta)=0$ has root $\bar{\lambda}$ with
${\rm{Re}}(\bar{\lambda})>0$. Then

$\left|{\bar{\lambda}(\bar{\lambda}+\eta(a+b))}\right|=\left|{{\eta^{2}}ab\bar{\lambda}\int\limits_{0}^{2}{{e^{-u\bar{\lambda}}}du}}\right|\leq{\eta^{2}}\left|{ab}\right||\bar{\lambda}|\left|{\int\limits_{0}^{2}{{e^{-u({\rm{Re}}(\bar{\lambda}))}}du}}\right|\leq
2{\eta^{2}}\left|{ab}\right||\bar{\lambda}|.$

Since $ab<0$ and $a+b+2\eta ab\geq 0$, we have

$|\bar{\lambda}(\bar{\lambda}+\eta(a+b))|\leq-2{\eta^{2}}ab|\bar{\lambda}|\leq\eta(a+b)|\bar{\lambda}|$

which is satisfied if $\bar{\lambda}=0$ (a contradiction) or

$|\bar{\lambda}+\eta(a+b)|\leq\eta(a+b).$

This implies that $\bar{\lambda}$ is in the disk of radius $\eta({a+b})$
centred at the point $-\eta({a+b})$ in the complex plane. Thus,
${\rm{Re}}(\bar{\lambda})<0$ or $\bar{\lambda}=0$. In both cases we arrive at
a contradiction. ∎

Finally, we show that (32) does not have pure imaginary roots for any value of
the parameters.

######  Proposition 3.6.

The characteristic equation (32) has no pure imaginary roots.

###### Proof.

Assume $\lambda=iy$ ($y>0$) is a root of (32). Separating the real and
imaginary parts, we obtain

$\displaystyle\eta^{2}ab-y^{2}$ $\displaystyle=\eta^{2}ab\cos(2y)$
$\displaystyle\eta(a+b)y$ $\displaystyle=-\eta^{2}ab\sin(2y)$

Squaring and adding these equations leads to

$y^{2}\left(y^{2}+\eta^{2}(a^{2}+b^{2})\right)=0.$

which has no real roots. Thus, there are no roots of the form $iy$. ∎

The distribution of roots in (32) is summarized in Figure 3.

Figure 3: The distribution of roots in (32) as discussed in Propositions
3.1$-$3.6.

Recall the structure of the phase-locked solutions (21) of the phase model
(20). From this we see that a phase-locked periodic solution of the original
model (17) corresponds to a line in the phase model (20), that is, when
$\psi^{*}$ and $\omega^{*}$ are solutions of (22), it follows that

$\left\\{{\begin{array}[]{*{20}{c}}{{\varphi_{1}}={\omega^{*}}t}&{\left({\bmod
2\pi}\right)}\\\ {{\varphi_{2}}={\omega^{*}}t+{\psi^{*}}}&{\left({\bmod
2\pi}\right)}\end{array}}\right.\Rightarrow{\varphi_{2}}={\varphi_{1}}+{\psi^{*}}\left({\bmod
2\pi}\right).$

From Proposition 3.4, we know that for any $\tau>0$,
$\Delta(\lambda;\eta(\tau))=0$ has a zero root. The simple zero root
corresponds to the motion along these lines. It corresponds to the Floquet
multiplier $1$ which is associated with the periodic solution of the original
model (17). Thus phase-locked solutions will be asymptotically stable if
$\lambda=0$ is a simple root of the characteristic equation (32) and all other
roots have negative real part.

######  Remark 3.1.

The solution $\psi^{*}\neq 0,\pi$ is asymptotically stable for values of $a,b$
such that $a>0$ and $b>0$ or $ab<0,\ a+b>0$ and $a+b+2\eta ab>0$. Since
$H^{\prime}$ is a $2\pi-$periodic function, the solutions $\psi^{*}$ and
$2\pi-\psi^{*}$ have the same stability.

######  Remark 3.2.

Since $H$ is a $2\pi-$periodic function,
$a=b=\frac{1}{\Omega}H^{\prime}(\psi^{*}-\omega^{*}\eta-\Omega\tau)$ in (30)
when $\psi^{*}=0,\pi$. Hence, the stability of solutions when $\psi^{*}=0,\pi$
is determined by the sign of $H^{\prime}(\psi^{*}-\omega^{*}\eta-\Omega\tau)$,
that is, the solution is asymptotically stable when
$H^{\prime}(\psi^{*}-\omega^{*}\eta-\Omega\tau)>0$ and unstable when
$H^{\prime}(\psi^{*}-\omega^{*}\eta-\Omega\tau)<0$.

###  3.1 Bifurcation

Suppose that $\Omega$ and $\epsilon$ are fixed, but $\tau$ may be varied. From
the discussion above, potential bifurcation points of the model (20) are
values $\tau=\tau^{*}$ where the characteristic equation for a particular
phase-locked solution, $\psi^{*},\omega^{*}$ has a double zero root. Let
$\eta^{*}=\epsilon\Omega\tau^{*}$. When $\psi^{*}=0$ or $\pi$ there are two
types of potential bifurcation points:

  * (1)

$\tau^{*}$ where $H^{\prime}(\psi^{*}-\omega^{*}\eta^{*}-\Omega\tau^{*})=0$
(see Remark 3.2);

  * (2)

$\tau^{*}$ where
$1+\eta^{*}\frac{1}{\Omega}H^{\prime}(\psi^{*}-\omega^{*}\eta^{*}-\Omega\tau^{*})=0$
(see Proposition 3.4).

For other values of $\psi^{*}$, Proposition 3.4 indicates there is a potential
bifurcation point at

  * (3)

$\tau^{*}$ where $\eta^{*}=-\frac{a+b}{2ab}$.

Note that it is impossible to find an explicit expression for the bifurcation
values because each of these conditions are implicit equations for $\tau^{*}$.

Now we consider what type of bifurcations may occur at these points. We do not
make a rigorous proof, which would require centre manifold and normal form
theory. However, we can make some plausible arguments based on the equations
for the equilibrium solutions. Recall that $(\psi^{*},\omega^{*})$ with
$\psi^{*}=0$ or $\pi$ defines a phase-locked solution at $\tau$ if
$F(\omega^{*},\psi^{*};\tau)=0$ where

$F(\omega,\psi^{*};\tau)=\omega-\frac{1}{\Omega}H(\psi^{*}-\omega\eta-\Omega\tau).$

Differentiating $F$ with respect to $\omega$ shows that the condition (2)
corresponds to $F_{\omega}(\omega^{*},\psi^{*};\tau^{*})=0$, that is,
$\omega^{*}$ is a double root of $F$ when $\tau=\tau^{*}$. Thus as $\tau$
varies near $\tau^{*}$ we may expect that there should be two roots of $F$
near $\omega^{*}$ or none 555More precisely, we expect this will occur if $F$
satisfies the further conditions
$F_{\tau}(\omega^{*},\psi^{*};\tau^{*})=-\Omega(1+\epsilon\omega^{*})/\eta^{*}\neq
0$ and
$F_{\omega\omega}(\omega^{*},\psi^{*};\tau^{*})=-(\eta^{*})^{2}H^{\prime\prime}(\psi^{*}-\omega^{*}\eta^{*}-\Omega\tau^{*})\neq
0$ [39].. Thus the bifurcation associated with condition (2) should be a
saddle-node bifurcation involving two different phase-locked solutions with
the same $\psi^{*}$. Note that this bifurcation is only physically relevant if
$\eta^{*}>0$, i.e.,
$H^{\prime}(\psi^{*}-\omega^{*}\eta^{*}-\Omega\tau^{*})<0$. Thus, from Remark
3.2, the associated solutions will be unstable. In a similar manner one can
show that condition (3) corresponds to $(\psi^{*},\omega^{*})$ at
$\tau=\tau^{*}$ being a point of tangency of the curves defined by equations
(22). Thus we expect it to correspond to a saddle-node bifurcation involving
two out-of-phase solutions with different $\psi^{*}$. The stability of these
solutions will depend on which case of Proposition 3.4 applies. Finally, we
consider phase-locked solutions near $\psi=0$. Expanding equations (25) and
the first of (22) in $\psi$ and keeping the two lowest order terms we have

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
2H^{\prime}(-\omega\eta-\Omega\tau)\psi+\frac{2}{3}H^{\prime\prime\prime}(-\omega\eta-\Omega\tau)\psi^{3}$
(36) $\displaystyle\omega$ $\displaystyle=$
$\displaystyle\frac{1}{\Omega}\left(H(-\omega\eta-\Omega\tau)+H^{\prime}(-\omega\eta-\Omega\tau)\psi\right).$
(37)

Thus we see that $\psi^{*}=0$,
$\omega^{*}=H(-\omega^{*}\eta-\Omega\tau)/\Omega$, is always a solution of
this system and if there is $\tau^{*}$ such that condition (1) is satisfied
and $H^{\prime\prime\prime}(-\omega^{*}\eta^{*}-\Omega\tau^{*})\neq 0$ then
this will be a triple root of the system. Thus we expect that condition (1)
with $\psi^{*}=0$ corresponds to a pitchfork bifurcation where two out-of-
phase solutions are created near $0$. Similarly condition (1) with
$\psi^{*}=\pi$ should correspond to a pitchfork bifurcation where two out-of-
phase solutions are created near $\pi$.

Note that the phase interaction function $H$ can be represented by Fourier
series expansion

$H(\phi)={a_{0}}+\sum\limits_{k=1}^{\infty}{\left[{{a_{k}}\cos(k\phi)+{b_{k}}\sin(k\phi)}\right]}.$

When the interaction function $H$ is represented by the first set of Fourier
modes

$H(\phi)=a_{0}+a_{1}\cos(\phi)+b_{1}\sin(\phi),$ (38)

the authors in [36] show that the out-of-phase solutions and pitchfork
bifurcation cannot occur in the phase model (20) with small time delay.
However, it may occur when the time delay is large. Indeed, when $H$ has the
form in (38), then it follows from (22) and (25) that

$\displaystyle\Omega{\omega^{*}}$
$\displaystyle={a_{0}}+A({\omega^{*}})\sin({\psi^{*}})+B({\omega^{*}})\cos({\psi^{*}}),$
(39) $\displaystyle 0$ $\displaystyle=2A({\omega^{*}})\sin({\psi^{*}})$ (40)

respectively, where

$\displaystyle A({\omega^{*}})$
$\displaystyle={b_{1}}\cos({\omega^{*}}\eta+\Omega\tau)+{a_{1}}\sin({\omega^{*}}\eta+\Omega\tau),$
$\displaystyle B({\omega^{*}})$
$\displaystyle={a_{1}}\cos({\omega^{*}}\eta+\Omega\tau)-{b_{1}}\sin({\omega^{*}}\eta+\Omega\tau).$

Thus, from $\sin({\psi^{*}})=0$, we have that $\psi^{*}=0,\pi$ with the
corresponding $\omega^{*}$ determined by

$\Omega{\omega^{*}}-{a_{0}}=\pm B({\omega^{*}}).$ (41)

where the $+$ corresponds to $\psi^{*}=0$ and the $-$ to $\psi^{*}=\pi$. Also,
from $A({\omega^{*}})=0$ we determine $\omega^{*}$ and the corresponding
$\psi^{*}$ is obtained from

$\cos({\psi^{*}})=\frac{{\Omega{\omega^{*}}-{a_{0}}}}{{B({\omega^{*}})}}.$
(42)

Consequently, we have the following cases

  * •

if $\left|{\Omega{\omega^{*}}-{a_{0}}}\right|<\left|{B({\omega^{*}})}\right|$,
then two out-of-phase solutions $\psi^{*}$ and $2\pi-\psi^{*}$ exist,

  * •

if $\left|{\Omega{\omega^{*}}-{a_{0}}}\right|=\left|{B({\omega^{*}})}\right|$,
then one solution exists ($\psi^{*}=0$ or $\psi^{*}=\pi$),

  * •

if $\left|{\Omega{\omega^{*}}-{a_{0}}}\right|>\left|{B({\omega^{*}})}\right|$,
then no solution satisfying (42) exists.

Note that $H^{\prime}({-\omega^{*}}\eta-\Omega\tau)=A(\omega^{*})$ and
$H^{\prime}({\pi-\omega^{*}}\eta-\Omega\tau)=-A(\omega^{*})$. Thus, the
solutions $0$ and $\pi$ change stability when $A(\omega^{*})=0$ where
$\omega^{*}$ satisfies (41). As $\tau$ varies, out-of-phase solutions will
disappear if $\frac{{\Omega{\omega^{*}}-{a_{0}}}}{{B({\omega^{*}})}}-1$
changes its sign from negative to positive. When
$\frac{{\Omega{\omega^{*}}-{a_{0}}}}{{B({\omega^{*}})}}=1$, then $\psi^{*}=0$.
Hence, a pitchfork bifurcation occurs at $\psi^{*}=0$. Similarly when
$\frac{{\Omega{\omega^{*}}-{a_{0}}}}{{B({\omega^{*}})}}=-1$ a pitchfork
bifurcation occurs at $\psi^{*}=\pi$.

###  3.2 The full model with small delay

When the time delay, $\tau$, in (17) is relatively small, in the sense that
$\Omega\tau=\mathcal{O}(1)$, it follows from the theory of averaging that the
time delay $\tau$ enters the interaction function $H$ in (20) as a phase shift
[24, 14, 27, 36]. In [36], the authors considered this case and consequently
the time delay $\eta$ in the phase model (20) was neglected, and hence, it
becomes

$\displaystyle\frac{d\varphi_{1}}{d{t}}$
$\displaystyle=\frac{1}{\Omega}H(\varphi_{2}(t)-\varphi_{1}(t)-\Omega\tau),$
(43) $\displaystyle\frac{d\varphi_{2}}{d{t}}$
$\displaystyle=\frac{1}{\Omega}H(\varphi_{1}(t)-\varphi_{2}(t)-\Omega\tau),$

Therefore, they were able to reduce (20) into a one dimensional ordinary
differential equation

$\frac{d\phi}{dt}=-2\epsilon[H(\phi-\Omega\tau)-H(-\phi-\Omega\tau)].$ (44)

where $\phi=\varphi_{2}-\varphi_{1}$. The existence of phase-locked solutions
of (44) was discussed in [36] without introducing the frequency deviation
$\omega$. Hence, the in-phase and anti-phase solutions were unique. Moreover,
the stability of the phase-locked solution $\phi^{*}$ in (44) was determined
by the sign of

${\widehat{H}}^{\prime}(\phi^{*}):=\overline{a}+\overline{b}$ (45)

where $\overline{a}=H^{\prime}(\phi^{*}-\Omega\tau)$ and
$\overline{b}=H^{\prime}(-\phi^{*}-\Omega\tau)$. If
${\widehat{H}}^{\prime}(\phi^{*})>0$ then $\phi^{*}$ is asymptotically stable
and if ${\widehat{H}}^{\prime}(\phi^{*})<0$ it is unstable. When
${\widehat{H}}^{\prime}(\phi^{*})=0$ the stability is not determined by the
linearization.

######  Remark 3.3.

In [36], due to the reduction of the two dimensional system (43) into a single
equation (44), the zero root was omitted in characteristic equation. Indeed,
the characteristic equation of (44) is
$\lambda+{\widehat{H}}^{\prime}(\phi^{*})=0$ while the characteristic equation
of (43) is

$\lambda({\lambda}+{\widehat{H}}^{\prime}(\phi^{*}))=0.$ (46)

It is clear that the latter characteristic equation always has a zero root.

Now we compare these results with what happens when $\tau$ is small, i.e.,
$\Omega\tau=\mathcal{O}(1)$, in our model (20). Recall that
$\eta=\epsilon\Omega\tau$ thus the assumption on $\tau$ implies that
$\eta=\mathcal{O}(\epsilon)$. Also, note that the phase difference $\phi^{*}$
of the phase locked solutions for the model (45) is the same as the phase
deviation difference $\psi^{*}$ for our model.

First consider the existence of phase-locked solutions. For our model we must
solve the equations (25) and one of (22) simultaneously for $\psi$ and
$\omega$. When $\eta=\mathcal{O}(\epsilon)$, however, to first order in
$\epsilon$ the $H$ function no longer depends on $\omega$. Thus phase-locked
solutions are determined by $\psi^{*}$ satisfying $H_{\tau}(\psi^{*})=0$, with
$\omega^{*}=\frac{1}{\Omega}H(\psi^{*}-\Omega\tau)$. This equation for
$\psi^{*}$ is the same as in [36]. In [36] they did not solve for $\omega^{*}$
as it was not needed to determine the phase-locked solutions or their
stability. It remains to consider the uniqueness of the in-phase and anti-
phase solutions. From equations (26) and (27), these solutions correspond to
frequency deviations $\omega^{*}$ satisfying $F(\omega^{*},\psi^{*})=0$ with
$\psi^{*}=0,\pi$, respectively. Since $H$ and $H^{\prime}$ are continuous and
$2\pi$ periodic they are bounded. Thus we see that
$\lim_{\omega\rightarrow\pm\infty}F(\omega)=\pm\infty$. Further, recalling
(28), since $\eta=\mathcal{O}(\epsilon)$, $F_{\omega}(\omega^{*},\psi^{*})>0$.
Thus for any $\tau$ sufficiently small, there will be a unique frequency
deviation $\omega^{*}$ for $\psi^{*}=0$ and for $\psi^{*}=\pi$. This is
consistent with the results in [36] which have only one in-phase and anti-
phase solution for each value of $\tau$.

Now consider the stability of the phase-locked solutions. Recall that the
stability for our model is summarized in Figure 3. When
$\eta=\mathcal{O}(\epsilon)$, ${\rm sgn}(a+b+ab\eta)\approx{\rm sgn}(a+b)$,
thus the conditions for stability/instability of phase-locked solutions of our
model reduce to the stability if $a+b>0$ and instability if $a+b<0$. Further
$a\approx\overline{a}$ and $b\approx\overline{b}$, thus the stability results
of our model reduce to those of [36] when $\Omega\tau=\mathcal{O}(1)$. The key
point is that, regardless of the size of $\tau$, the countable infinity of
complex roots of the characteristic equation (32) all have negative real part.
Thus the stability of the phase-locked solutions is determined by finitely
many real roots, and it is possible for an ordinary differential equation to
accurately reflect this stability.

In Section 4.3, we will show numerically that our model with
$\Omega\tau=\mathcal{O}(1)$ fully recovers [36, Figure 4b] and [36, Figure
5b].

##  4 Application to Morris-Lecar oscillators with diffusive
coupling

In this section we apply the results from the previous sections to a network
of dimensionless Morris-Lecar oscillators with time delayed diffusive
coupling, see e.g., [40, 41]. This model is given by

$\displaystyle{v^{\prime}_{i}}$
$\displaystyle={I_{app}}-{g_{Ca}}{m_{\infty}}({v_{i}})({v_{i}}-{v_{Ca}})-{g_{K}}{w_{i}}({v_{i}}-{v_{K}})-{g_{L}}({v_{i}}-{v_{L}})-\epsilon({v_{j}}(t-\tau)-{v_{i}}(t)),$
(47) $\displaystyle{w^{\prime}_{i}}$
$\displaystyle=\varphi\lambda({v_{i}})({w_{\infty}}({v_{i}})-{w_{i}}),$

for $i,j=1,2$ such that $i\neq j$, where

$\displaystyle m_{\infty}(v)$
$\displaystyle=\frac{1}{2}\left(1+\tanh\left(\left(v-\nu_{1}\right)/\nu_{2}\right)\right),$
$\displaystyle w_{\infty}(v)$
$\displaystyle=\frac{1}{2}\left(1+\tanh\left(\left(v-\nu_{3}\right)/\nu_{4}\right)\right),$
$\displaystyle\lambda(v)$
$\displaystyle=\cosh\left(\left(v-\nu_{3}\right)/\left(2\nu_{4}\right)\right).$

Using the parameter set I$\backslash$II from [36, Table 1], when there is no
coupling in the network each oscillator has a unique exponentially
asymptotically stable limit cycle with period $T=23.87\backslash 13.81$
corresponding to frequency $\Omega=0.2632\backslash 0.455$. The normalized
system, such that the frequency is $1$, corresponding to (47) is

$\displaystyle{v^{\prime}_{i}}$
$\displaystyle=\frac{1}{\Omega}({I_{app}}-{g_{Ca}}{m_{\infty}}({v_{i}})({v_{i}}-{v_{Ca}})-{g_{K}}{w_{i}}({v_{i}}-{v_{K}})-{g_{L}}({v_{i}}-{v_{L}}))-\frac{\epsilon}{\Omega}({v_{j}}(t-\Omega\tau)-{v_{i}}(t)),$
(48) $\displaystyle{w^{\prime}_{i}}$
$\displaystyle=\frac{1}{\Omega}(\varphi\lambda({v_{i}})({w_{\infty}}({v_{i}})-{w_{i}})),$

$i=1,2$. Note that this is in the form (17) with
${\bf{X}}_{i}(t)=(v_{i}(t),w_{i}(t))^{T}$ and the function
$\mathbf{G}:\mathbb{R}^{2}\times\mathbb{R}^{2}\to\mathbb{R}^{2}$ is given by
${\bf{G}}=(G_{1},G_{2})$ where
${{G}}_{1}({\bf{X}}_{1}(t),{\bf{X}}_{2}(t))=\frac{1}{\Omega}(v_{2}(t−\Omega\tau)−v_{1}(t))$
and ${{G}}_{2}({\bf{X}}_{1}(t),{\bf{X}}_{2}(t))=0$. Then, the phase model
interaction function $H$ is given by (19).

For each parameter set, the authors in [36] solved (19) numerically and
calculated the approximation of the phase model interaction function $H$ by
the first five terms of its Fourier series. These are given by

$\displaystyle H_{I}(\phi)$
$\displaystyle=2.915252-2.684797\cos(\phi)-0.3278022\cos(2\phi)$ (49)
$\displaystyle~{}~{}~{}+0.05596774\cos(3\phi)+0.0351635\cos(4\phi)+4.908449\sin(\phi)$
$\displaystyle~{}~{}~{}-0.7020183\sin(2\phi)-0.09934668\sin(3\phi)-0.01104474\sin(4\phi),$
$\displaystyle H_{II}(\phi)$
$\displaystyle=0.6271561-0.5209326\cos(\phi)-0.08538575\cos(2\phi)$
$\displaystyle~{}~{}~{}-0.005648281\cos(3\phi)-0.0002642404\cos(4\phi)+1.595618\sin(\phi)$
$\displaystyle~{}~{}~{}-0.04727176\sin(2\phi)-0.00301241\sin(3\phi)-0.002760313\sin(4\phi)$

corresponding to the parameter sets I and II, respectively, see [36, Table 2].
Note that the two parameter sets represent limit cycles which are created by
different bifurcations as the input current $I_{app}$ is varied. For parameter
set I the limit cycle is created in a saddle-node on an invariant circle
bifurcation, while for parameter set II the limit cycle is created in a
supercitical Hopf bifurcation. The chosen parameter values have $I_{app}$
slightly larger than the bifurcation values.

In [36] the authors studied how small epsilon needed to be for the phase model
to faithfully represented the behaviour of the full system (48), in the case
of small delay. They found that for parameter set I $\epsilon$ could be as
large as $0.05$ while for parameter set II epsilon should not exceed $0.001$.
Therefore, in the rest of this section, we take $\epsilon=0.05$ with parameter
set I and $\epsilon=0.001$ when we use parameter set II. Consequently, we
choose $\tau\geq 75.988$ for parameter set I and $\tau\geq 2197.8$ for
parameter set II so that $\epsilon\Omega\tau=\mathcal{O}(1)$. Moreover, we
compare our results with the results in [36] when the time delay, $\tau$, in
(17) is relatively small.

###  4.1 In-phase and anti-phase solutions

To find $\omega^{*}$ corresponding to the in-phase and anti-phase solutions,
$\psi^{*}=0,\pi$, we solve (26) and (27) with $H$ given by either $H_{I}$ or
$H_{II}$ from (49). Note that these equations can only be solved numerically
due to the complicated form of $H_{I}$ and $H_{II}$. For particular values of
$\tau$, we represent these solutions graphically in Figure 4 as the
intersection points of the line $y=\omega$ and the curve
$y=H(-\omega\eta-\Omega\tau)/\Omega$. In (26), the slope of the right hand
side at any $\omega$ is $\ell_{\tau}=-\epsilon\tau
H^{\prime}(-\omega\eta-\Omega\tau)$. Then, by applying the stability condition
in Remark 3.2, we see that the in-phase solution is stable when the line
$y=\omega$ intersects the curve of the function
$H(-\omega\eta-\Omega\tau)/\Omega$ at a point where it has negative slope,
while it is unstable when the intersection is at a point with positive slope,
see Figure 4. When the line $y=\omega$ alternates from intersecting the curve
of $H(-\omega\eta-\Omega\tau)/\Omega$ at a point with positive slope to
intersecting it at a point with negative slope, the solutions $\omega^{*}$
alternate between stable and unstable, see Figure 4. For fixed $\Omega$ and
$\epsilon$, as $\tau$ increases the curve
$H(-\eta\omega-\Omega\tau)=H(-\Omega\tau(1+\epsilon\omega))$ compresses
horizontally causing the creation and destruction of intersection points. For
specific values $\tau=\tau^{*}_{1}>0$, an intersection point will occur at the
point where the function $H$ has slope one, i.e., the curve $y=H$ will be
tangent to the line $y=\omega$ at these values of $\tau$, see Figure 5(b).
Near such points, i.e., for $\tau$ slightly bigger or smaller, there exist two
consecutive intersection points both of which are unstable, see Figure 5(c).
Then, as $\tau$ changes further to $\tau_{2}^{*}$, one unstable point quickly
passes through the point where $H$ has zero slope and becomes stable, see
Figure 4(b). The values $\tau^{*}_{1}$ correspond to the saddle-node
bifurcations of in-phase and anti-phase solutions discussed in Section 3.1. We
will discuss the points $\tau_{2}^{*}$ later. In Figure 7, we plot
$\omega^{*}$ corresponding to $\psi^{*}=0,\pi$ for various values of the time
delay $\tau$, showing the many co-existing solutions which can occur and the
transitions of the solutions as $\tau$ varies. These solutions were found by
implementing the algorithm from [42] in Wolfram Mathematica to find all the
solutions of (26) or (27).

(a) $H_{I}(-\omega\eta-\Omega\tau)/\Omega$ and $y=\omega$. $\tau=90$.

(b) $H_{II}(-\omega\eta-\Omega\tau)/\Omega$ and $y=\omega$. $\tau=2236$.

Figure 4: Graphical representation of the solutions to (26) with fixed $\tau$.
The circles  $\CIRCLE$/ $\Circle$ represent stable/unstable solutions.

(a) $\tau=2234.4$.

(b) $\tau=2234.78$.

(c) $\tau=2235$.

Figure 5: Graphical solutions of $H_{II}(-\omega\eta-\Omega\tau)/\Omega$ and
$y=\omega$ with fixed $\tau$. The circles  $\CIRCLE$/$\Circle$ represent
stable/unstable solutions.

To compare prediction of the phase model (20) and solutions of the full model
(47), we solve (48) numerically with parameter sets I and II with various
values of $\tau$ and different initial conditions. The initial conditions are
of the form

$\displaystyle\left(v_{1}(t),w_{1}(t),v_{2}(t),w_{2}(t)\right)^{T}=\left(v_{10},w_{10},v_{20},w_{20}\right)^{T}\quad
t\in[-\tau\Omega,0].$ (50)

Figure 6 shows time series of $v_{i}$ in (48) with different initial
conditions. We notice the coexistence of in-phase solutions with different
frequencies when $\tau=110$ with parameter set I. The numerical solutions are
obtained by using Wolfram Mathematica. We use the command NDSolve to solve the
full model numerically.

Figure 6: The coexistence of in-phase solutions of (48) with different
frequencies when $\tau=110$ with parameter sets I. We take different initial
conditions: $(1.53422,-4.42364,1.58103,-4.12258)^{T}$ for red/orange curves
and $(-2.22807,3.626,-2.28885,-0.972632)^{T}$ for blue/green curves.

When $\epsilon=0$, each uncoupled equation in (48) has $2\pi-$periodic
solution, that is, the frequency of each oscillator is unity. Consequently,
when $\epsilon\neq 0$ and equation (48) has a phase-locked solution, the phase
of the first oscillator is $\theta_{1}(t)=t+\omega^{*}\epsilon t$ and that of
the second oscillator is $\theta_{2}(t)=t+\omega^{*}\epsilon t+\psi^{*}$ where
$\omega^{*}$ is the frequency deviation and $\psi^{*}$ is the phase shift.
Thus, the frequency of each oscillator is $1+\omega\epsilon^{*}$, and the
period $\mathcal{T}$ is approximately

$\mathcal{T}=\frac{2\pi}{1+\omega^{*}\epsilon}.$

From the numerical solution of (48) for a stable phase-locked solution, we can
calculate the period $\mathcal{T}$ of the oscillators and determine the
approximate frequency deviation from

$\displaystyle\omega^{*}\approx\frac{1}{\epsilon}\left(\frac{2\pi}{\mathcal{T}}-1\right).$
(51)

Figure 7 shows the coexistence of stable in-phase and anti-phase periodic
solutions and demonstrates that the approximation of $\omega^{*}$ from (51) is
close to a stable solution of the phase model. The values of $\omega^{*}$ with
the normalized error

${{\rm
E_{N}}}=\frac{(\omega^{*}{\rm{~{}in~{}the~{}phase~{}model}})-(\omega^{*}{\rm{~{}in~{}the~{}full~{}model}})}{\omega^{*}{\rm{~{}in~{}the~{}full~{}model}}}$
(52)

are shown in Tables 1 and 2. Note that the quantity ${{\rm E_{N}}}$ is the
normalized error with respect to the size of $\omega^{*}$ in the full model.
Except for a few cases, the phase model gives a very accurate prediction of
the values of $\omega^{*}$. The phase model predicted stable phase-locked
solutions that we did not find numerically, however, it is possible that
further exploration with different initial conditions might find them.

(a) Parameter set I. $\psi^{*}=0$.

(b) Parameter set I. $\psi^{*}=\pi$

(c) Parameter set II. $\psi^{*}=0$

(d) Parameter set II. $\psi^{*}=\pi$

Figure 7:  The circles $\CIRCLE$/$\Circle$ represent stable/unstable solutions
to the phase model (20) corresponding to (48) and $\boldsymbol{\times}$
represents the calculated $\omega^{*}$ for each stable phase-locked periodic
solution found by numerical integration of the full model (48) with parameter
sets I and II. The insets are in-phase and anti-phase periodic solutions of
(48). The initial conditions for all simulations was of the form (50). For the
insets the values of $\left(v_{10},w_{10},v_{20},w_{20}\right)^{T}$ are as
follows. (a) $(1.53422,-4.42364,1.58103,-4.12258)^{T}$ (b)
$(1.57882,4.1827,2.78262,0.358165)^{T}$ (c)
$(1.892,-0.296437,-1.05518,1.09985)^{T}$ (d)
$(-1.72448,-1.46442,4.4848,1.31822)^{T}$

$\psi^{*}=0$
---
| $\tau=90$ |  | $\tau=110$ |  | $\tau=130$
Phase Model | Full Model | ${\rm E_{N}}$ |  | Phase Model | Full Model | ${\rm E_{N}}$ | Phase Model | Full Model | ${\rm E_{N}}$
$\omega^{*}$ | $1.1446$ | $0.988971$ | $0.1574$ |  | $1.61521$ | $1.44374$ | $0.1188$ | $-1.55286$ | $−1.31666$ | $-0.1794$
$6.17015$ | $5.68396$ | $0.0855$ |  | $9.95867$ | $9.67259$ | $0.0296$ | $5.48587$ | $5.16611$ | $0.0619$
| $\tau=150$ |  | $\tau=170$ |  | $\tau=190$
Phase Model | Full Model | ${\rm E_{N}}$ |  | Phase Model | Full Model | ${\rm E_{N}}$ | Phase Model | Full Model | ${\rm E_{N}}$
$\omega^{*}$ | $-0.860951$ | $−0.766851$ | $-0.1227$ |  | $2.38636$ | $2.23525$ | $0.0676$ | $0.101849$ | $0.092197$ | $0.1047$
$2.19521$ | $2.03459$ | $0.0789$ |  | $5.11589$ | $4.87837$ | $0.0487$ | $7.44601$ | $7.16109$ | $0.0398$
$\psi^{*}=\pi$
| $\tau=90$ |  | $\tau=110$ |  | $\tau=130$
Phase Model | Full Model | ${\rm E_{N}}$ |  | Phase Model | Full Model | ${\rm E_{N}}$ | Phase Model | Full Model | ${\rm E_{N}}$
$\omega^{*}$ | $-1.32864$ | $-1.08976$ | $-0.2192$ |  | $3.68511$ | $3.39257$ | $0.0862$ | $0.193408$ | $0.169843$ | $0.1387$
$8.70987$ | $8.244$ | $0.0565$ |  | $7.86$ | $7.36178$ | $0.0677$ | $7.26467$ | $6.86785$ | $0.0578$
| $\tau=150$ |  | $\tau=170$ |  | $\tau=190$
Phase Model | Full Model | ${\rm E_{N}}$ |  | Phase Model | Full Model | ${\rm E_{N}}$ | Phase Model | Full Model | ${\rm E_{N}}$
$\omega^{*}$ | $0.663926$ | $0.600602$ | $0.1054$ |  | $1.02812$ | $0.946137$ | $0.0867$ | $3.76194$ | $3.58466$ | $0.0495$
$9.92786$ | $9.41424$ | $0.0546$ |  | $3.74932$ | $3.55154$ | $0.0557$ | $8.67729$ | $8.34016$ | $0.0404$

Table 1: Comparison of $\omega^{*}$ between the phase model prediction and the
full model (48) when $\psi^{*}=0,\pi$ with parameter set I. The quantity ${\rm
E_{N}}$ is defined in (52).

$\psi^{*}=0$
---
| $\tau=2200$ |  | $\tau=2500$
Phase Model | Full Model | ${\rm E_{N}}$ |  | Phase Model | Full Model | ${\rm E_{N}}$
$\omega^{*}$ | $-1.418638$ | $-1.10245$ | $-0.2868$ |  | $-0.159684$ | $-0.223359$ | $0.2851$
| $\tau=2800$ |  | $\tau=3100$
Phase Model | Full Model | ${\rm E_{N}}$ |  | Phase Model | Full Model | ${\rm E_{N}}$
$\omega^{*}$ | $0.9584587$ | $0.720141$ | $0.3309$ |  | $1.914776$ | $1.880915$ | $0.018$
$\psi^{*}=\pi$
| $\tau=2200$ |  | $\tau=2500$
Phase Model | Full Model | ${\rm E_{N}}$ |  | Phase Model | Full Model | ${\rm E_{N}}$
$\omega^{*}$ | $0.913646$ | $0.647528$ | $0.411$ |  | $2.06122$ | $1.58187$ | $0.303$
| $\tau=2800$ |  | $\tau=3100$
Phase Model | Full Model | ${\rm E_{N}}$ |  | Phase Model | Full Model | ${\rm E_{N}}$
$\omega^{*}$ | $-1.011297$ | $-0.880915$ | $-0.148$ |  | $0.0496789$ | $-0.03801$ | $-2.307$

Table 2: Comparison of $\omega^{*}$ between the phase model prediction and the
full model (48) when $\psi^{*}=0,\pi$ with parameter set II. The quantity The
quantity ${\rm E_{N}}$ is defined in (52).

###  4.2 Out-of-phase solutions

To find phase-locked solutions other than the in-phase and anti-phase
solutions, we fix $\tau$ and solve

$\displaystyle\omega^{*}$
$\displaystyle=\frac{1}{\Omega}H_{II}(\psi^{*}-\omega^{*}\eta-\Omega\tau),$
(53) $\displaystyle\omega^{*}$
$\displaystyle=\frac{1}{\Omega}H_{II}(-\psi^{*}-\omega^{*}\eta-\Omega\tau)$

for $\omega^{*}$ and $\psi^{*}$. Figure 8 shows all solutions to (53) when
$\tau=100\backslash 2205$ with the parameter set I$\backslash$II. As seen for
the existence of in-phase and anti-phase solutions in Section 4.1, the number
of phase-locked solutions with the parameter set I is bigger that II. For the
purpose of clarity in the bifurcation figures, we consider the parameter set
II in this section.

In Figure 8(b), we observe that there are four non-trivial phase-locked
solutions: $\psi^{*}_{1}=1.85996$ and $\psi^{*}_{2}=2.13981$ in $(0,\pi)$; and
$\psi^{*}_{3}=2\pi-\psi^{*}_{1}=4.42323$ and
$\psi^{*}_{4}=2\pi-\psi^{*}_{2}=4.14338$ in $(\pi,2\pi)$. Moreover, we have
$\omega^{*}_{1}=\omega^{*}_{3}=0.14125$ and
$\omega^{*}_{2}=\omega^{*}_{4}=0.14125$ where $\omega^{*}_{i}$ is the
corresponding frequency deviation to $\psi^{*}_{i}$, $i=1,2,3,4$. This agrees
with Proposition 2.1.

(a) Parameter set I with $\tau=100$.

(b) Parameter set II with $\tau=2205$.

Figure 8: Contour plot of the equations in (53) to show the graphical
solutions of (53) with fixed $\tau$.

In Figure 9a, we plot all solutions of system (53) in $\tau\psi-$plane and
mark the stability using the criteria in Section 3. Note that since this
representation suppresses $\omega^{*}$, the multiple in-phase or anti-phase
solutions which occur for particular values of $\tau$ in Figure 7 are
superimposed. As $\tau$ varies, we observe that a stable solution
corresponding to $\psi^{*}=0,\pi$ always exists with the appearance of an
unstable solution in disjoint intervals of $\tau$, while all the out-of-phase
solutions are unstable. More precisely, for the in-phase solution, as $\tau$
increases, we notice that an unstable solution disappears at $\tau\approx
2203$, exists between $\tau\approx 2207.5$ and $\tau\approx 2217$, and
reappears at $\tau\approx 2221$. The same behaviour occurs for the anti-phase
solution at different values of $\tau$. Near the appearance and disappearance
of these unstable solutions the unstable out-of-phase solutions appear and
disappear. As we observe in Figure 8, there are multiple solutions
$(\omega^{*},\psi^{*})$ of (53) when $\tau$ is fixed. To study the creation
and destruction of solutions further, we take particular values for $\tau$ and
show all solutions in the blue rectangles from Figure 9a in the
$\omega\psi-$plane, see Figures 9b$-$9i. We now see that there are pitchfork
bifurcations where a stable in-phase or anti-phase solution becomes unstable
as two unstable out-of-phase solutions merge together, see Figures 9b$-$9c and
9f$-$9g. This correspond to the values $\tau_{2}^{*}$ discussed above.
Moreover, there are saddle-node bifurcations where two unstable in-phase or
anti-phase solutions collide then vanish, see Figures 9d$-$9e and 9h$-$9i.
This corresponds to the value $\tau_{1}^{*}$ discussed above. For other
parameter values, we observe the opposite sequence of bifurcations: two
unstable in-phase or anti-phase solutions are created by a saddle-node
bifurcation after which one gets stabilized by a pitchfork bifurcation
involving two unstable out-of-phase solutions. All the bifurcations are as
predicted for the general model in Section 3.1. We did not observe any saddle-
node bifurcations of out-of-phase solutions for this parameter set.

To help understand these bifurcations, we plot solutions in the
$\tau\omega-$plane and the solutions near $\psi=\pi$ in the
$\tau\omega\psi-$space in Figures 10(a)$-$10(b), respectively. Considering the
case $\psi^{*}=\pi$, we observe that:

  * •

the pitchfork bifurcation occurs when two unstable out-of-phase solutions
merge together with one stable anti-phase solution $\CIRCLE$ to produce one
unstable anti-phase solution $\blacksquare$,

  * •

the saddle-node bifurcation occurs when the created unstable anti-phase
solution $\blacksquare$ in the above collides with another unstable anti-phase
$\blacksquare$ and both vanish.

Figure 9: The solutions of phase model (20) corresponding to (48) in the blue
rectangles in Figure 10 in $\omega\psi-$plane. The circles $\CIRCLE$/$\Circle$
represent stable/unstable solutions of (53).

(a) $\omega$ vs $\tau$

(b) $\tau$, $\omega$ and $\psi$

Figure 10: Numerical bifurcation diagram with respect to $\tau\in(2200,2225)$
for the solutions of the phase model (20) corresponding to the Morris-Lecar
model (48) with parameter set II. The circles $\CIRCLE$/❂ represent
stable/unstable in-phase solutions, $\blacktriangle$/$\blacksquare$ represents
stable/unstable anti-phase solutions, and $\times$ represents unstable out-of-
phase solutions of (53).

###  4.3 Small delay

In this subsection, we consider small time delay, in the sense that,
$\Omega\tau=\mathcal{O}(1)$ with respect to the small parameter $\epsilon$,
and compare the results with [36] where the authors studied this case using
the parameter set II. In [36], the authors studied the dynamics of the phase
model corresponding to the full model (17) without introducing the frequency
deviation in their analysis because the time delay $\eta$ was neglected in the
phase model when $\Omega\tau=\mathcal{O}(1)$. We have stated some results from
[36] in Section 3.2.

As in the previous section we solve (26) and (27) to find $\omega^{*}$ for the
in-phase and anti-phase solutions and (53) to find $(\psi^{*},\omega^{*})$ for
the out-of-phase solutions. We choose $\tau\in(0,15)$, which is similar to the
range chosen by [36]. In contrast with the results of the last section, here
we observe that for $\psi^{*}=0,\pi$ there is a unique solution $\omega^{*}$
for each $\tau$ in the range we considered. This agrees with the prediction of
the phase model in Section 3.2. We describe our results in more detail below.

(a) $\tau\in(0,15)$

(b) $\tau\in(2.7,2.84)$

(c) $\tau\in(9.59,9.73)$

Figure 11:  Numerical bifurcation diagram with respect to $\tau\in(0,15)$ for
the phase model (20) corresponding to the Morris-Lecar model (48) with
parameter set II. The circles $\CIRCLE$/$\Circle$ represent stable/unstable
solutions in the phase model (53).

In Figure 11, we plot the in-phase and anti-phase solutions as $\tau$ varies
in $(0,15)$ in the $\tau\psi-$plane. We note that there is similar behaviour
in Figure 11(a) and [36, Figure 4b]. The in-phase and anti-phase solutions
change stability as $\tau$ increases and their stabilities appear to be the
opposite of each other. To examine the behaviour near changes of stability, in
Figures 11(b)$-$11(c) we show the bifurcation diagrams zoomed close to the two
switching points. We see that the transition from stable in-phase solution to
stable anti-phase solution involves two pitchfork bifurcations and one saddle-
node bifurcation of out-of-phase solutions, which agrees with [36]. Figure 13
shows this behaviour when the solutions are plotted in the $\tau\omega-$plane.
Furthermore, we observe in Figures 11(b)$-$11(c) that there are small
intervals of $\tau$ where bistability occurs. Figure 12 shows the coexistence
of stable anti-phase and out-of-phase solutions.

Figure 12: The coexistence of stable anti-phase (red) and out-of-phase (blue)
solutions of (48) when $\tau=9.661$ with parameter sets II. We take different
initial conditions: $(0.664192,0.204054,5.58914,0.762568)^{T}$ for red curve
and $(-0.883364,-0.200879,-0.686477,-0.989329)^{T}$ for blue curve.

######  Remark 4.1.

The results in this section are consistent with the results in [27], which
indicate that a phase model where the time delay enters as a phase shift is
accurate when $\tau$ is small in the full model (17) in the sense that
$\Omega\tau=\mathcal{O}(1)$ with respect to $\epsilon$ for $0<\epsilon\ll 1$.

(a) $\tau\in(0,15)$

(b) $\tau\in(2.7,2.84)$

(c) $\tau\in(9.59,9.73)$

Figure 13: Numerical bifurcation diagram with respect to $\tau$ for the phase
model (20) corresponding to the Morris-Lecar model (48) with parameter set II.
The circles $\CIRCLE$/$\Circle$ represent stable/unstable solutions in the
phase model (53).

##  5 Conclusions

In this paper, we studied the phase-locking dynamics of a system of two weakly
connected oscillators with time-delayed interaction. By applying the theory of
weakly coupled oscillators, we transformed the system into a phase model with
an explicit delay in the argument of the phases. We showed that the system
always has phase-locked solutions corresponding to in-phase (synchronous, $0$
phase difference) and anti-phase (phase difference of half the period)
solutions. Further, we showed for small delay ($\Omega\tau=\mathcal{O}(1)$)
the in-phase and anti-phase solutions are unique, but for large delay multiple
solutions of each type may exist, corresponding to different frequencies.
Finally, we showed that phase-locked solutions with any other phase
differences (out-of-phase solutions) are also possible. Since the phase model
is an infinite-dimensional system of delay differential equations, the
linearized system about the phase-locked solutions has a countable infinity of
eigenvalues. Through the stability analysis for our model, we discussed the
distribution of the eigenvalues on the complex plane to provide stability
conditions for the in-phase, anti-phase and out-of-phase solutions. We found
that the zero eigenvalue always exists for any choice of parameters and
functions which corresponds to the motion along the phase-locked solutions. We
showed that the only way in which bifurcations can occur is through the
existence of (additional) zero eigenvalues and argued that the following
bifurcations may occur: saddle-node bifurcations of two in-phase solutions
with different frequencies, saddle-node bifurcations of two anti-phase
solutions with different frequencies, saddle-node bifurcations of two
different out-of-phase solutions, pitchfork bifurcations where two out-of-
phase solutions arise from an in-phase or anti-phase solution. We showed that
the saddle-node bifurcations of in-phase and anti-phase solutions only involve
unstable solutions.

Our results on in-phase and anti-phase solutions agree with those in [23, 24],
which study the phase model (3), with $n=2$ and $H(\cdot)=\sin(\cdot)$. We
note that they emphasized the need for large coupling-strength for multiple
in-phase/anti-phase solutions to exist, however, we show that it is possible
with weak coupling and sufficiently large delays. They do not study out-of-
phase solutions as these are not possible in their model due to the
restriction on $H$. As can be seen in the literature [8, 9, 11, 12], in order
for phase models derived from biophysical oscillator models to adequately
capture the dynamics of the full model, the function $H$ generally must
include multiple Fourier modes. In [36] it was shown that out-of-phase
solutions and pitchfork bifurcations cannot occur in a phase model with small
delay if only the first Fourier modes are included in $H$. However, when the
time delay is large, we showed that both out-of-phase solutions and pitchfork
bifurcations can occur in the phase model with only the first Fourier modes of
$H$. In general, in the case of large time delay, the bifurcation structure
may change if some modes are dropped. If the coefficients of the modes that
are dropped are small, then the bifurcation structure wouldn’t change much.
The bifurcation points may just move around. If the coefficients of the modes
dropped are big enough then there could be large changes in the bifurcation
structure.

When the delay is small ($\Omega\tau=\mathcal{O}(1)$), Campbell and
Kobelevskiy studied the system

$\displaystyle\frac{d\theta_{1}}{d{t}}$ $\displaystyle=\Omega+\epsilon
H(\theta_{2}(t)-\theta_{1}(t)-\Omega\tau),$ (54)
$\displaystyle\frac{d\theta_{2}}{d{t}}$ $\displaystyle=\Omega+\epsilon
H(\theta_{1}(t)-\theta_{2}(t)-\Omega\tau),$

and proved that in-phase and anti-phase solutions are stable when
$H^{\prime}(\phi^{*}-\Omega\tau)>0$, $\phi^{*}\in\\{0,\pi\\}$ in [36]. On the
other hand, when the time delay is large $\epsilon\Omega\tau=\mathcal{O}(1)$,
we proved that these solutions are stable whenever
$H^{\prime}(\phi^{*}-\omega^{*}\epsilon\Omega\tau-\Omega\tau)>0$ where
$\omega^{*}$ is the corresponding frequency deviation. It is clear that the
stability condition in the first case is independent of the coupling strength
parameter and the frequency deviation. Indeed, under the assumption
$\theta_{1}(t)=\Omega+\omega t$ and $\theta_{2}(t)=\Omega+\omega t+\phi^{*}$
(see (21)), the terms of the frequency deviation $\omega$ will cancel out
inside the function $H$ in (54). In fact, in [36], the authors reduce (54)
into a single ordinary differential equation and study the dynamics of the
model without introducing the frequency deviation. Due to the explicit delay
in the phase model, we couldn’t reduce the model into a single equation. For
the out-of-phase solutions $\phi^{*}\notin\\{0,\pi\\}$, the stability
condition $H^{\prime}(\phi^{*}-\Omega\tau)>0$ is still valid when the delay is
small. While for the large delay the stability becomes more complicated since
the explicit delay is an additional parameter that needs to be considered in
the phase model.

As an example we considered two Morris-Lecar oscillators with delayed,
diffusive coupling. We adopted the parameter values from [36] to compare the
results when the time delay is small. We studied the existence and stability
of the phase-locked solutions, and explored the bifurcations in the phase
model by using a four mode trunction of the Fourier series for the interaction
function and compared these results with numerical simulations of the full
model. When the time delay $\tau$ is large, we found:

  * •

There exist more than one frequency deviation $\omega$ corresponding to the
in-phase and anti-phase solutions, i.e., co-existence of multiple stable and
unstable solutions;

  * •

All out-of-phase solutions are unstable;

  * •

Both the pitchfork and saddle-node bifurcations of in-phase and anti-phase
solutions occur.

When the time delay is small, we observed:

  * •

Unique solution in each phase-locked solution category (in-phase, anti-phase
and out-of-phase).

  * •

The occurrence of saddle-node bifurcations of out-of-phase solutions and
pitchfork bifurcations of in-phase and anti-phase solutions.

Our results agree with [36] when the time delay is small and are consistent
with the results in [27], that the explicit time delay can be neglected in the
phase model when $\tau$ is small.

A special type of phase-locked solutions, so-called symmetric cluster
solutions, can appear in a network of $n$ identical oscillators, see e.g.,
[17, 43],

$\frac{d\mathbf{X}_{i}}{dt}={\mathbf{F}}\left(\mathbf{X}_{i}(t)\right)+\epsilon\sum_{j=1}^{n}a_{ij}{\mathbf{G}}\left(\mathbf{X}_{i}(t),\mathbf{X}_{j}\left(t-\tau\right)\right),\quad
i=1,\ldots,n,\quad\mathbf{X}_{i}\in\mathbb{R}^{m}.$ (55)

In these solutions, also called travelling wave solutions, oscillators in the
same cluster are synchronized while those in different clusters have non-zero
phase-difference. In [17], Campbell and Wang determined conditions for
existence and stability of symmetric cluster solutions in (55) when $\tau$ is
small and the coupling matrix is circulant. Stability conditions for cluster
solutions in networks with small distance dependent delays and random, nearest
neighbour coupling have been formulated by several authors (see [24, 37] and
references therein). When the time delay is large, Earl and Strogatz provided
the stability condition for the in-phase solution ($\theta_{i}(t)=\Omega t$,
i.e., one cluster solution), see [18]. For future research, it would be
interesting to study the existence and stability of symmetric cluster
solutions in (55) with large time delay.

## Acknowledgments

The authors would like to thank the anonymous referees for their careful
reading and helpful suggestions.

##  Appendix A Phase reduction

Assume that the system (10) admits an exponentially asymptotically stable
periodic orbit with natural frequency $\Omega$ when $\epsilon=0$. It follows
from the time rescaling $\rho\to\Omega\rho$ that the natural frequency of the
periodic orbit becomes $1$ and (10) can be written as

$\frac{{d{{\mathbf{X}}_{i}}}}{{d\rho}}=\frac{1}{\Omega}{\mathbf{F}}({{\mathbf{X}}_{i}}(\rho))+\frac{\epsilon}{\Omega}\sum\limits_{j=1}^{n}{{K_{ij}}}{\mathbf{G}}({{\mathbf{X}}_{i}}(\rho),{{\mathbf{X}}_{j}}(\rho-\Omega\tau)),\quad
i=1,\ldots n,\ {{\mathbf{X}}_{i}}\in{\mathbb{R}^{m}}$ (A.1)

Consequently, there exists a normally hyperbolic invariant manifold
$M=\gamma\times\cdots\times\gamma$ of system (A.1) when $\epsilon=0$, where
$\gamma$ is an exponentially orbitally stable $2\pi-$periodic solution of

$\frac{d\mathbf{X}_{i}}{d\rho}=\frac{1}{\Omega}{\mathbf{F}}({\mathbf{X}_{i}}(\rho))\quad
i=1,\ldots,n.$ (A.2)

Hence, the solution of the $i^{th}$ equation of (A.1) in an $\epsilon$
neighborhood of $M$ can be written as

${{\mathbf{X}}_{i}}(\rho)=\gamma\left({\rho+{\varphi_{i}}(t)}\right)+\epsilon{P_{i}}\left({\rho,\varphi_{1}(t),\ldots,\varphi_{n}(t),\epsilon}\right)$
(A.3)

where the term $\epsilon{P_{i}}$ is a smooth vector function which denotes the
deviation from the manifold $M$ in the normal plane.

Recall that $t=\epsilon\rho$ and let $\eta:=\epsilon\Omega\tau$, then the
substitution of (A.3) in (A.1) gives

$\displaystyle\frac{{d{{\mathbf{X}}_{i}}}}{{d\rho}}$
$\displaystyle=\frac{1}{\Omega}{\mathbf{F}}\left[{\gamma\left({\rho+{\varphi_{i}}(t)}\right)+\epsilon
P_{i}\left({\rho,\varphi_{1}(t),\ldots,\varphi_{n}(t),\epsilon}\right)}\right]$
(A.4)
$\displaystyle~{}+\frac{\epsilon}{\Omega}\sum_{j=1}^{n}K_{ij}{\mathbf{G}}\left[{\gamma\left({\rho+{\varphi_{i}}(t)}\right)+\epsilon
P_{i}\left({\rho,\varphi_{1}(t),\ldots,\varphi_{n}(t),\epsilon}\right),}\right.$
$\displaystyle~{}\qquad\qquad\left.{\gamma\left({\rho-\Omega\tau+{\varphi_{j}}(t-\eta)}\right)+\epsilon
P_{j}\left({\rho,{\varphi_{1}}(t-\eta),\ldots,{\varphi_{n}}(t-\eta),\epsilon}\right)}\right].$

Due to the infinite differentiability of ${\mathbf{F}}$ and ${\mathbf{G}}$, it
follows from (A.4) that

$\displaystyle\frac{{d{{\mathbf{X}}_{i}}}}{{d\rho}}=\frac{1}{\Omega}{\mathbf{F}}\left[{\gamma\left({\rho+{\varphi_{i}}(t)}\right)}\right]+\frac{\epsilon}{\Omega}D{\mathbf{F}}\left[{\gamma\left({\rho+{\varphi_{i}}(t)}\right)}\right]P_{i}\left({\rho,\varphi_{1}(t),\ldots,\varphi_{n}(t),\epsilon}\right)$
(A.5)
$\displaystyle+\frac{\epsilon}{\Omega}\sum_{j=1}^{n}K_{ij}{\mathbf{G}}\left[{\gamma\left({\rho+{\varphi_{i}}(t)}\right),\gamma\left({\rho-\Omega\tau+{\varphi_{j}}(t-\eta)}\right)}\right]+\mathcal{O}\left({{\epsilon^{2}}}\right)$

where $D{\mathbf{F}}$ is the Jacobian matrix of ${\mathbf{F}}$.

Now, we differentiate ${\mathbf{X}}_{i}$ in (A.3) with respect to $\rho$ to
have

$\frac{{d{{\mathbf{X}}_{i}}}}{{d\rho}}=\gamma^{\prime}\left({\rho+{\varphi_{i}}(t)}\right)\left({1+\epsilon\frac{{d{\varphi_{i}}}}{{dt}}}\right)+\epsilon\frac{{\partial{P_{i}}\left({\rho,\varphi_{1}(t),\ldots,\varphi_{n}(t),\epsilon}\right)}}{{\partial\rho}}+\mathcal{O}\left({{\epsilon^{2}}}\right).$
(A.6)

Note that

${\gamma^{\prime}}\left({\rho+{\varphi_{i}}(t)}\right)=\frac{1}{\Omega}{\mathbf{F}}\left[{\gamma\left({\rho+{\varphi_{i}}(t)}\right)}\right].$
(A.7)

Thus, from (A.5) and (A.6), we obtain

$\displaystyle{\mathbf{F}}\left[{\gamma\left({\rho+{\varphi_{i}}(t)}\right)}\right]\frac{{d{\varphi_{i}}(t)}}{{dt}}$
$\displaystyle+\Omega\frac{{\partial
y_{i}\left({\rho,\varphi_{1}(t),\ldots,\varphi_{n}(t)}\right)}}{{\partial\rho}}=D{\mathbf{F}}\left[{\gamma\left({\rho+{\varphi_{i}}(t)}\right)}\right]y\left({\rho,\varphi_{1}(t),\ldots,\varphi_{n}(t)}\right)$
$\displaystyle{+\sum_{j=1}^{n}K_{ij}{\mathbf{G}}\left[{\gamma\left({\rho+{\varphi_{i}}(t)}\right),\gamma\left({\rho-\Omega\tau+{\varphi_{j}}(t-\eta)}\right)}\right]+\mathcal{O}(\epsilon)}$
(A.8)

where
$y_{i}\left({\rho,\varphi_{1}(t),\ldots,\varphi_{n}(t)}\right):=P_{i}\left({\rho,\varphi_{1}(t),\ldots,\varphi_{n}(t),0}\right)+\mathcal{O}(\epsilon)$
because $P_{i}$ is smooth function of $\epsilon$. Consequently, since
$t=\epsilon\rho$, we replace $\frac{\partial y_{i}}{\partial\rho}$ by
$\frac{dy_{i}}{d\rho}$ in (A). Hence, we can write (A) as:

$\frac{{dy_{i}}}{{d\rho}}=A\left({\rho,{\varphi_{i}}}\right)y+{b_{i}}\left({\rho,{\varphi_{1}},\ldots,{\varphi_{n}}}\right)+\mathcal{O}(\epsilon)$
(A.9)

where $\varphi_{i}$ is $\varphi_{i}(t)$,

$A\left({\rho,{\varphi_{i}}}\right)=\frac{1}{\Omega}{\mathbf{F}}\left[{\gamma\left({\rho+{\varphi_{i}}(t)}\right)}\right]$

and

${b_{i}}\left({\rho,{\varphi_{i}},\ldots,{\varphi_{n}}}\right)=\frac{1}{\Omega}\left[{\left({\sum\limits_{j=1}^{n}{{K_{ij}}}{\bf{G}}\left[{\gamma\left({\rho+{\varphi_{i}}(t)}\right),\gamma\left({\rho-\Omega\tau+{\varphi_{j}}(t-\eta)}\right)}\right]}\right)-{\bf{F}}\left[{\gamma\left({\rho+{\varphi_{i}}(t)}\right)}\right]\frac{{d{\varphi_{i}}(t)}}{{dt}}}\right].$

Since ${{\varphi_{i}}(t)}$ and ${{\varphi_{i}}(t-\eta)}$ in $b$ do not depend
directly on $\rho$, we have a linear non-homogeneous system for $y_{i}$, where
both the matrix $A$ and the vector $b_{i}$ are $2\pi$ periodic in $\rho$.

To study existence and uniqueness of solutions to (A.9), we consider the
adjoint linear homogeneous system

$\frac{{d{Q_{i}{\left({\rho,{\varphi_{i}}}\right)}}}}{{d\rho}}=-A{\left({\rho,{\varphi_{i}}}\right)^{T}}{Q_{i}}{\left({\rho,{\varphi_{i}}}\right)}$
(A.10)

with the normalization condition:

$\frac{1}{{2\pi}}\int\limits_{0}^{2\pi}{Q_{i}^{T}\left({\rho,{\varphi_{i}}}\right)}{\mathbf{F}}\left[{\gamma\left({\rho+{\varphi_{i}}}\right)}\right]d\rho=1.$
(A.11)

Since the limit cycle $\gamma$ is exponentially orbitally stable, the
homogeneous ($b_{i}\equiv 0$) linear system of the form (A.9) the adjoint
system (A.10) both have $1$ as a simple Floquet multiplier, and all the other
multipliers lie inside the unit circle. Thus, system (A.10)-(A.11) has a
unique nontrivial periodic solution
${q_{i}\left({\rho,{\varphi_{i}}}\right)}$.

Now, by the Fredholm alternative, the linear non-homogeneous system (A.9) has
a unique periodic solution $Y_{i}$ if and only if the following orthogonality
condition holds:

$\left\langle{{q_{i}},{b_{i}}}\right\rangle+\mathcal{O}(\epsilon)=\frac{1}{{2\pi}}\int\limits_{0}^{2\pi}{q_{i}^{T}\left({\rho,{\varphi_{i}}}\right)}{b_{i}}\left({\rho,{\varphi_{i}},\ldots,{\varphi_{n}}}\right)d\rho+\mathcal{O}(\epsilon)=0.$
(A.12)

Assume that $q_{i}(\rho,0)$ is found. Hence,
$q_{i}(\rho,\varphi_{i})=q_{i}(\rho+\varphi_{i},0)$ because
$A\left({\rho,{\varphi_{i}}}\right)=\frac{1}{\Omega}D{\mathbf{F}}\left[{\gamma\left({\rho+{\varphi_{i}}(t)}\right)}\right]=A\left({\rho+{\varphi_{i}}(t)}\right)$.
Thus, when we substitute $b_{i}$ in (A.12), we obtain the following:

$\displaystyle\frac{1}{{2\pi\Omega}}\int\limits_{0}^{2\pi}{q_{i}^{T}\left({\rho+{\varphi_{i}},0}\right)}\left({\sum\limits_{j=1}^{n}{{K_{ij}}{\bf{G}}\left[{\gamma\left({\rho+{\varphi_{i}}(t)}\right),\gamma\left({\rho-\Omega\tau+{\varphi_{j}}(t-\eta)}\right)}\right]}}\right)d\rho+{\cal
O}(\epsilon)$ $\displaystyle\qquad\hskip
170.71652pt=\frac{1}{{2\pi}}\int\limits_{0}^{2\pi}{q_{i}^{T}\left({\rho+{\varphi_{i}},0}\right)}D{\mathbf{F}}\left[{\gamma\left({\rho+{\varphi_{i}}(t)}\right)}\right]\frac{{d{\varphi_{i}}(t)}}{{dt}}d\rho.$

Since $\frac{{d{\varphi_{i}}(t)}}{{dt}}$ is treated as a parameter and is
independent of $\rho$, it follows from the normalization condition (A.11) that

$\frac{{d{\varphi_{i}}(t)}}{{dt}}=\frac{1}{{2\pi\Omega}}\int\limits_{0}^{2\pi}{q_{i}^{T}\left({\rho+{\varphi_{i}},0}\right)}\left({\sum\limits_{j=1}^{n}{{K_{ij}}{\bf{G}}\left[{\gamma\left({\rho+{\varphi_{i}}(t)}\right),\gamma\left({\rho-\Omega\tau+{\varphi_{j}}(t-\eta)}\right)}\right]}}\right)d\rho+{\cal
O}(\epsilon)$

Letting $s=\rho+\varphi_{i}$ leads to

$\frac{{d{\varphi_{i}}(t)}}{{dt}}=\frac{1}{{2\pi\Omega}}\sum\limits_{j=1}^{n}{{K_{ij}}}\left({\int\limits_{0}^{2\pi}{q_{i}^{T}\left({s,0}\right)}{\bf{G}}\left[{\gamma\left(s\right),\gamma\left({s-\Omega\tau+{\varphi_{j}}(t-\eta)-{\varphi_{i}}(t)}\right)}\right]ds}\right)+{\cal
O}(\epsilon)$

Define

$H\left({{\varphi_{j}}(t-\eta)-{\varphi_{i}}(t)-\Omega\tau}\right)=\frac{1}{{2\pi}}\int\limits_{0}^{2\pi}{q_{i}^{T}\left({s,0}\right)}{\mathbf{G}}\left[{\gamma\left(s\right),\gamma\left({s-\Omega\tau+{\varphi_{j}}(t-\eta)-{\varphi_{i}}(t)}\right)}\right]ds.$

Thus, we have system (11) with

$\frac{{d{\varphi_{i}}(t)}}{{dt}}=\frac{1}{{\Omega}}\sum\limits_{j=1}^{n}{{K_{ij}}}H\left({{\varphi_{j}}(t-\eta)-{\varphi_{i}}(t)-\Omega\tau}\right)+\mathcal{O}(\epsilon).$
(A.13)

Recall that $\eta=\epsilon\Omega\tau$. Hence, when $\Omega\tau=\mathcal{O}(1)$
with respect to $\epsilon$, we have

$\varphi_{i}(t-\eta)=\varphi_{i}(t-\epsilon\Omega\tau)=\varphi_{i}(t)+\mathcal{O}(\epsilon),\qquad
i=1,\ldots,n.$

Consequently, the Taylor series expansion for $h$ with respect to $\epsilon$
gives

$\displaystyle
H\left({{\varphi_{j}}(t-\epsilon\Omega\tau)-{\varphi_{i}}(t)-\Omega\tau}\right)$
$\displaystyle=H\left({{\varphi_{j}}(t)-{\varphi_{i}}(t)-\Omega\tau+\mathcal{O}(\epsilon)}\right)$
$\displaystyle=H\left({{\varphi_{j}}(t)-{\varphi_{i}}(t)-\Omega\tau}\right)+\mathcal{O}(\epsilon),$

that is, no delay appears in the argument of the phases. Hence, (A.13)
becomes:

$\frac{{d{\varphi_{i}}(t)}}{{dt}}=\frac{1}{{\Omega}}\sum\limits_{j=1}^{n}{{K_{ij}}}H\left({{\varphi_{j}}(t)-{\varphi_{i}}(t)-\Omega\tau}\right)+\mathcal{O}(\epsilon).$
(A.14)

## References

  * [1] Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence.  New York: Springer-Verlag, 1984.
  * [2] R. E. Mirollo and S. H. Strogatz, “Synchronization of pulse-coupled biological oscillators,” SIAM J. Appl. Math., vol. 50, no. 6, pp. 1645–1662, 1990\.
  * [3] H. G. Winful and S. S. Wang, “Stability of phase locking in coupled semiconductor laser arrays,” Appl. Phys. Lett., vol. 53, no. 20, pp. 1894–1896, 1988.
  * [4] S. S. Wang and H. G. Winful, “Dynamics of phase-locked semiconductor laser arrays,” Appl. Phys. Lett., vol. 52, no. 21, pp. 1774–1776, 1988.
  * [5] F. Dörfler, M. Chertkov, and F. Bullo, “Synchronization in complex oscillator networks and smart grids,” Proceedings of the National Academy of Sciences, vol. 110, no. 6, pp. 2005–2010, 2013.
  * [6] N. Kopell and G. B. Ermentrout, “Coupled oscillators and the design of central pattern generators,” Math. Biosci., vol. 90, no. 1-2, pp. 87–109, 1988\.
  * [7] D. Hansel, G. Mato, and C. Meunier, “Phase dynamics for weakly coupled Hodgkin-Huxley neurons,” Europhys. Lett., vol. 23, no. 5, pp. 367–372, 1993.
  * [8] S. M. Crook, G. B. Ermentrout, M. C. Vanier, and J. M. Bower, “The role of axonal delay in the synchronization of networks of coupled cortical oscillators,” Journal of computational neuroscience, vol. 4, no. 2, pp. 161–172, 1997.
  * [9] Y. Park and B. Ermentrout, “Weakly coupled oscillators in a slowly varying world,” Journal of computational neuroscience, vol. 40, no. 3, pp. 269–281, 2016.
  * [10] A. Takamatsu, T. Fujii, and I. Endo, “Time delay effect in a living coupled oscillator system with plasmodium of physarum polycephalum,” Phys. Rev. E, vol. 85, no. 9, pp. 2026–2029, 2000.
  * [11] E. Wall, F. Guichard, and A. R. Humphries, “Synchronization in ecological systems by weak dispersal coupling with time delay,” Theoretical ecology, vol. 6, no. 4, pp. 405–418, 2013.
  * [12] Y. X. Zhang, F. Lutscher, and F. Guichard, “How robust is dispersal-induced spatial synchrony?,” Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 25, no. 3, p. 036402, 2015.
  * [13] F. Dörfler and F. Bullo, “Synchronization in complex networks of phase oscillators: a survey,” Automatica, vol. 50, pp. 1539–1564, 2014.
  * [14] F. C. Hoppensteadt and E. M. Izhikevich, Weakly Connected Neural Networks, vol. 126.  New York: Springer-Verlag, 1997.
  * [15] M. Porter and J. P. Gleeson, Dynamical Systems on Networks: A Tutorial.  Springer, 2016.
  * [16] M. A. Schwemmer and T. J. Lewis, “The theory of weakly coupled oscillators,” in Phase Response Curves in Neuroscience (N. W. Schultheiss, A. A. Prinz, and R. J. Butera, eds.), pp. 3–31, New York: Springer, 2012.
  * [17] S. A. Campbell and Z. Wang, “Phase models and clustering in networks of oscillators with delayed coupling,” Physica D, vol. 363, pp. 44–55, 2018\.
  * [18] M. G. Earl and S. H. Strogatz, “Synchronization in oscillator networks with delayed coupling: A stability criterion,” Phys. Rev. E, vol. 67, no. 3, p. 036204, 2003.
  * [19] S. Kim, S. H. Park, and C. S. Ryu, “Multistability in coupled oscillator systems with time delay,” Phys. Rev. Lett., vol. 79, no. 15, pp. 2911–2914, 1997.
  * [20] T. Luzyanina, “Synchronization in an oscillator neural network model with time-delayed coupling,” Network: Computation in Neural Systems, vol. 6, pp. 43–59, 1995.
  * [21] E. Niebur, H. Schuster, and D. Kammen, “Collective frequencies and metastability in networks of limit-cycle oscillators with time delay,” Phys. Rev. Lett., vol. 67, pp. 2753–2756, 1991.
  * [22] M. K. Yeung and S. H. Strogatz, “Time delay in the Kuramoto model of coupled oscillators,” Phys. Rev. Lett., vol. 82, no. 3, pp. 648–651, 1999.
  * [23] H. Schuster and P. Wagner, “Mutual entrainment of two limit cycle oscillators with time delayed coupling,” Prog. Theor. Phys., vol. 82, no. 5, pp. 939–945, 1989.
  * [24] B. Ermentrout and T. W. Ko, “Delays and weakly coupled neuronal oscillators,” Philos. Trans. R. Soc. A-Math. Phys. Eng. Sci., vol. 367, no. 1891, pp. 1097–1115, 2009.
  * [25] M. Brede and A. C. Kalloniatis, “Frustration tuning and perfect phase synchronization in the Kuramoto-Sakaguchi model,” Phys. Rev. E, vol. 93, no. 6, p. 062315, 2016.
  * [26] H. Sakaguchi and Y. Kuramoto, “A soluble active rotater model showing phase transitions via mutual entertainment,” Prog. Theor. Phys., vol. 76, no. 3, pp. 576–581, 1986.
  * [27] E. M. Izhikevich, “Phase models with explicit time delays,” Phys. Rev. E, vol. 58, no. 1, pp. 905–908, 1998.
  * [28] B. Ermentrout, “An introduction to neural oscillators,” in Neural modeling and neural networks (F. Ventriglia, ed.), pp. 79–110, Oxford, UK: Pergamon Press, 1994.
  * [29] G. Ermentrout and D. Terman, Mathematical Foundations of Neuroscience.  New York, NY: Springer, 2010.
  * [30] N. Kopell and G. Ermentrout, “Mechanisms of phase-locking and frequency control in pairs of coupled neural oscillators,” in Handbook of Dynamical Systems, vol 2: Toward Applications (B. Fiedler, ed.), pp. 3–54, Amsterdam: Elsevier, 2002.
  * [31] R. F. Galán, “The phase oscillator approximation in neuroscience: an analytical framework to study coherent activity in neural networks,” in Coordinated Activity in the Brain, pp. 65–89, Springer, 2009.
  * [32] T. Zahid and F. Skinner, “Predicting synchronous and asynchronous network groupings of hippocampal interneurons coupled with dendritic gap junctions.,” Brain Research, vol. 1262, pp. 115–129, 2009.
  * [33] B. Pietras and A. Daffertshofer, “Network dynamics of coupled oscillators and phase reduction techniques,” Phys. Rep., 2019.
  * [34] H. Nakao, “Phase reduction approach to synchronisation of nonlinear oscillators,” Contemp. Phys., vol. 57, no. 2, pp. 188–214, 2016.
  * [35] P. Ashwin, S. Coombes, and R. Nicks, “Mathematical frameworks for oscillatory network dynamics in neuroscience,” J. Math. Neurosci., vol. 6, no. 1, p. 2, 2016.
  * [36] S. A. Campbell and I. Kobelevskiy, “Phase models and oscillators with time delayed coupling,” Discret. Contin. Dyn. Syst. Ser. A, vol. 8, pp. 2653–2673, 2012.
  * [37] T.-W. Ko, S.-O. Jeong, and H.-T. Moon, “Wave formation by time delays in randomly coupled oscillators,” Phys. Rev. E, vol. 69, no. 5, p. 056106, 2004.
  * [38] A. Pikovsky and M. Rosenblum, “Synchronization,” Scholarpedia, vol. 2, no. 12, p. 1459, 2007.
  * [39] Y. A. Kuznetsov, Elements of Applied Bifurcation Theory.  New York: Springer, 1998.
  * [40] A. Prasad, S. K. Dana, R. Karnatak, J. Kurths, B. Blasius, and R. Ramaswamy, “Universal occurrence of the phase-flip bifurcation in time-delay coupled systems,” Chaos, vol. 18, no. 2, p. 023111, 2008.
  * [41] N. Burić and D. Todorović, “Dynamics of Fitzhugh-Nagumo excitable systems with delayed coupling,” Phys. Rev. E, vol. 67, no. 6, p. 066222, 2003.
  * [42] S. K. Rahimian, F. Jalali, J. Seader, and R. E. White, “A new homotopy for seeking all real roots of a nonlinear equation,” Computers & chemical engineering, vol. 35, no. 3, pp. 403–411, 2011.
  * [43] K. Okuda, “Variety and generality of clustering in globally coupled oscillators,” Physica D, vol. 63, pp. 424–436, 1993.