Isochronous Dynamics in Pulse Coupled Oscillator Networks with Delay

TL;DR

This paper investigates how delay-induced discontinuities in pulse-coupled oscillator networks lead to the formation of isochronous regions filled with periodic orbits, providing analytical and numerical insights into their properties.

Contribution

It introduces the concept of isochronous regions in pulse-coupled oscillator networks with delay and offers a detailed analytical and numerical analysis of their existence and characteristics.

Findings

Identification of isochronous regions with same-period orbits

Analytical proof of isochronous region existence

Numerical characterization of multiple isochronous regions

Abstract

We consider a network of identical pulse-coupled oscillators with delay and all-to-all coupling. We demonstrate that the discontinuous nature of the dynamics induces the appearance of isochronous regions---subsets of the phase space filled with periodic orbits having the same period. For fixed values of the network parameters each such isochronous region corresponds to a subset of initial states on an appropriate surface of section with non-zero dimension such that all periodic orbits in this set have qualitatively similar dynamical behaviour. We analytically and numerically study in detail such an isochronous region, give a proof of its existence, and describe its properties. We further describe other isochronous regions that appear in the system.

Figures

b](#S2.F2.sf2). Note that the function $H$ in Eq. ([3](#S2.E3)) has the property

b](#S3.F3.sf2). We denote the corresponding isochronous region by IR4.

a](#S4.F4.sf1) and the detailed dynamics is given in [Table 1](#S4.T1). The set $\Omega_{\varepsilon,\tau}$ also gives rise to periodic orbits with smaller minimal orbit period than $T=3\tau$. In particular, there is a line in $\Omega_{\varepsilon,\tau}$ given by $(\sigma_{1},\sigma_{2},\sigma_{3})=(\tau/2,\sigma_{2},\tau/2+\sigma_{2})$ for which all points give rise to period $T=3\tau/2$ orbits ($T_{P}=2$), see [Fig. 4b](#S4.F4.sf2). One point along this line, having $\sigma_{2}=\tau/4$ gives rise to a period $T=3\tau/4$ orbit ($T_{P}=1$), see [Fig. 4c](#S4.F4.sf3).

a](#S4.F5.sf1) for $(\varepsilon,\tau)=(0.58,0.58)$. One can check that $\Omega_{\varepsilon,\tau}$ and $\mathcal{A}_{\varepsilon,\tau}$ have non-empty interior for $(\varepsilon,\tau)=(0.58,0.58)$ and that each point in $\widetilde{\Omega}_{\varepsilon,\tau}$ is the initial condition of a periodic orbit in IR4 of [Table 1](#S4.T1) with period $T=3\tau$ and Poincaré period $T_{P}=4$ . Therefore, $\widetilde{\Omega}_{\varepsilon,\tau}$ is a periodic plateau.

a](#S4.F7.sf1) compares the analytically obtained $\mathcal{A}_{\varepsilon,\tau}$ for $(\varepsilon,\tau)=(0.58,0.58)$ to the numerical results discussed in [subsection III.1](#S3.SS1). We note that the numerically obtained orbits cover only part of $\mathcal{A}_{\varepsilon,\tau}$. This can be explained by the fact that the space of initial conditions that we scanned in our numerical experiments does not include the periodic orbits in IR4. Some of the orbits in IR4 are periodic attractors for our numerical initial conditions but others are not accessible. This effect is much more pronounced

a](#S5.F10.sf1).

a](#S5.F12.sf1).

Tables

Table 1. Table 1: Dynamics for a periodic orbit in the isochronous region IR4. The periodic orbits in this isochronous region have Poincaré period T P = 4 subscript 𝑇 𝑃 4 T_{P}=4 and period T = 3 ​ τ 𝑇 3 𝜏 T=3\tau in the full phase space. In the time interval ( − τ , 0 ) 𝜏 0 (-\tau,0) , the only useful information is the firing moments, so we use ‘ − - ’ in the table to represent the “useless” information and ‘ F 𝐹 F ’ to represent an oscillator fired at the given moment.

time	$O_1[\left\{{\sigma }_1\right\},{\theta }_1$]	$O_2[\left\{{\sigma }_2\right\},{\theta }_2$]	$O_3[\left\{{\sigma }_3\right\},{\theta }_3$]
$-{\sigma }_3$	$[-,-]$	$[-,-]$	$[-,F]$
$-{\sigma }_1$	$[-,F]$	$[-,-]$	$[-,-]$
$-{\sigma }_2$	$[-,-]$	$[-,F]$	$[-,-]$
0	$[{\sigma }_1,H\left({\sigma }_1\right)]$	$[{\sigma }_2,{\sigma }_2]$	$[{\sigma }_3,0]$
$\tau -{\sigma }_3$	$[\tau -{\sigma }_3+{\sigma }_1,0]$	$[\tau -{\sigma }_3+{\sigma }_2,H\left(\tau -{\sigma }_3+{\sigma }_2\right)]$	$[\tau -{\sigma }_3,\tau -{\sigma }_3]$
$\tau -{\sigma }_1$	$[{\sigma }_3-{\sigma }_1,{\sigma }_3-{\sigma }_1]$	$[\tau -{\sigma }_1+{\sigma }_2,0]$	$[\tau -{\sigma }_1,H\left(\tau -{\sigma }_1\right)]$
$\tau -{\sigma }_2$	$[{\sigma }_3-{\sigma }_2,H\left({\sigma }_3-{\sigma }_2\right)]$	$[{\sigma }_1-{\sigma }_2,{\sigma }_1-{\sigma }_2]$	$[\tau -{\sigma }_2,0]$
$\tau$	$[{\sigma }_3,0]$	$[{\sigma }_1,H\left({\sigma }_1\right)]$	$[{\sigma }_2,{\sigma }_2]$
$2\tau -{\sigma }_3$	$[\tau -{\sigma }_3,\tau -{\sigma }_3]$	$[\tau -{\sigma }_3+{\sigma }_1,0]$	$[\tau -{\sigma }_3+{\sigma }_2,H\left(\tau -{\sigma }_3+{\sigma }_2\right)]$
$2\tau -{\sigma }_1$	$[\tau -{\sigma }_1,H\left(\tau -{\sigma }_1\right)]$	$[{\sigma }_3-{\sigma }_1,{\sigma }_3-{\sigma }_1]$	$[\tau -{\sigma }_1+{\sigma }_2,0]$
$2\tau -{\sigma }_2$	$[\tau -{\sigma }_2,0]$	$[{\sigma }_3-{\sigma }_2,H\left({\sigma }_3-{\sigma }_2\right)]$	$[{\sigma }_1-{\sigma }_2,{\sigma }_1-{\sigma }_2]$
$2\tau$	$[{\sigma }_2,{\sigma }_2]$	$[{\sigma }_3,0]$	$[{\sigma }_1,H\left({\sigma }_1\right)]$
$3\tau -{\sigma }_3$	$[\tau -{\sigma }_3+{\sigma }_2,H\left(\tau -{\sigma }_3+{\sigma }_2\right)]$	$[\tau -{\sigma }_3,\tau -{\sigma }_3]$	$[\tau -{\sigma }_3+{\sigma }_1,0]$
$3\tau -{\sigma }_1$	$[\tau -{\sigma }_1+{\sigma }_2,0]$	$[\tau -{\sigma }_1,H\left(\tau -{\sigma }_1\right)]$	$[{\sigma }_3-{\sigma }_1,{\sigma }_3-{\sigma }_1]$
$3\tau -{\sigma }_2$	$[{\sigma }_1-{\sigma }_2,{\sigma }_1-{\sigma }_2]$	$[\tau -{\sigma }_2,0]$	$[{\sigma }_3-{\sigma }_2,H\left({\sigma }_3-{\sigma }_2\right)]$
$3\tau$	$[{\sigma }_1,H\left({\sigma }_1\right)]$	$[{\sigma }_2,{\sigma }_2]$	$[{\sigma }_3,0]$

Equations

$${\dot{x}}_i(t)$$

$$x_i(t^{+})$$

$$x_i(t)$$

$$x_i=f({\theta }_i),$$

$$\dot{x}=F(x)=f^{\prime }(f^{-1}(x))=\frac{1}{(f^{-1}{)}^{\prime }(x)},$$

$$F(x):=F_b(x)=\frac{e^b-1}{b}{e}^{-bx},b>0,$$

$$f(\theta ):=f_b(\theta )=\frac{1}{b}\ln \left(1+(e^b-1)\theta \right).$$

$${\dot{\theta }}_i(t)$$

$${\theta }_i(t^{+})$$

$${\theta }_i(t)$$

$$H(\theta ,\delta )=f^{-1}(f(\theta )+\delta )=e^{b\delta }\theta +\frac{e^{b\delta }-1}{e^b-1},$$

$$V(\theta ,\delta )=H(\theta ,\delta )-\theta =(e^{b\delta }-1)\theta +\frac{e^{b\delta }-1}{e^b-1},$$

$$H(H(\theta ,\delta ),{\delta }^{\prime })=H(\theta ,\delta +{\delta }^{\prime }).$$

$$H(H(\theta ,m\hat{\epsilon }),m^{\prime }\hat{\epsilon })=H(\theta ,(m+m^{\prime })\hat{\epsilon }).$$

$$H(\theta ,m\hat{\epsilon })=H_m(\theta )\text{and}H(\theta ,\hat{\epsilon })=H_1(\theta )=H(\theta ).$$

$$\theta :[-\tau ,0]\to {\mathbb{T}}^n:t\mapsto ({\theta }_1(t),\dots ,{\theta }_N(t)).$$

$$\cdots <-{\sigma }_i^{(3)}<-{\sigma }_i^{(2)}<-{\sigma }_i^{(1)}\le 0.$$

$$\displaystyle\Phi=\bigl{\{}\{\sigma_{i}^{(j)}\}_{j},\theta_{i}\bigr{\}}_{i=1,\dots,N}.$$

$$\displaystyle\Phi=\bigl{\{}\{\sigma_{i}^{(j)}\}_{j},\theta_{i}\bigr{\}}_{i=1,\dots,N}\mapsto\{\theta_{i}\}_{i=1,\dots,N-1}.$$

$${\Sigma }_i=\{{\sigma }_i^{(j)}\}=\{\begin{array}{ll}\{{\theta }_i\},& \text{if θi≤τ}\\ \varnothing ,& \text{if θi>τ}.\end{array}$$

$$\begin{array}{rl}G& :({\sigma }_1,{\sigma }_2,{\sigma }_3;{\theta }_1,{\theta }_2)\mapsto \\ & ({\sigma }_3-{\sigma }_2,{\sigma }_1-{\sigma }_2,\tau -{\sigma }_2;H({\sigma }_3-{\sigma }_2),{\sigma }_1-{\sigma }_2).\end{array}$$

$$[P,(1,2),\tau -{\sigma }_3],[P,(2,3),\tau -{\sigma }_1],[P,(1,3),\tau -{\sigma }_2],$$

$$[P,(1,2),\tau ];[F,1,1-{\theta }_1],[F,2,1-{\theta }_2],[F,3,1].$$

$$\stackrel{1}{⟶}$$

$$[P,(1,2),{\sigma }_3];[F,1,1-{\theta }_1-\tau +{\sigma }_3],$$

$$[F,2,1-{\theta }_2-\tau +{\sigma }_3],[F,3,1-\tau +{\sigma }_3]$$

$$\stackrel{2}{⟶}$$

$$[F,1,0],[F,2,1-H({\theta }_2+\tau -{\sigma }_3)],[F,3,1-\tau +{\sigma }_3].$$

$$\stackrel{3}{⟶}$$

$$[P,(1,2),{\sigma }_3],[P,(2,3),\tau ];[F,2,1-H({\theta }_2+\tau -{\sigma }_3)],$$

$$[F,3,1-\tau +{\sigma }_3],[F,1,1]$$

$$\stackrel{4}{⟶}$$

$$[P,(2,3),\tau +{\sigma }_1-{\sigma }_3];$$

$$[F,2,1-H({\theta }_2+\tau -{\sigma }_3)+{\sigma }_1-{\sigma }_3],$$

$$[F,3,1-\tau +{\sigma }_1],[F,1,1+{\sigma }_1-{\sigma }_3]$$

$$\stackrel{5}{⟶}$$

$$[P,(2,3),\tau +{\sigma }_1-{\sigma }_3];[F,2,0],$$

$$[F,3,1-H(\tau -{\sigma }_1)],[F,1,1+{\sigma }_1-{\sigma }_3]$$

$$\stackrel{6}{⟶}$$

$$[P,(2,3),\tau +{\sigma }_1-{\sigma }_3],[P,(1,3),\tau ];$$

Full text

pulse coupled oscillator networks, isochronous dynamics, synchronization

Pulse coupled oscillator networks are a key model for the study of synchronization in a wide variety of systems, ranging from fireflies to wireless communication systems. Moreover, despite their simplicity, they manifest dynamical behavior that does not typically appear in smooth finite-dimensional dynamical systems. We report on the existence of isochronous dynamics in pulse coupled oscillator networks with delay: for suitable values of the parameters there exist open sets of initial conditions giving periodic orbits with the same period. This, previously unknown, behavior of pulse coupled oscillator networks with delay provides a deeper understanding of their dynamics and how they can reach synchronization.

I Introduction

Pulse coupled oscillator networks

Pulse coupled oscillator networks (PCONs) have been used to model interactions in networks where each node affects other nodes in a discontinuous way. Two such examples are the synchronization related to the function of the heart (Peskin, 1975) and the synchronization of fireflies (Mirollo and Strogatz, 1990). There is now an extensive literature on the dynamics of pulse coupled oscillator networks focusing on synchronization and the stability of synchronized states.

Concerning syncronization, after the seminal work Mirollo and Strogatz (1990) who considered excitatory coupling with no delay, Ernst, Pawelzik, and Geisel (1995, 1998) showed the importance of delayed and inhibitory coupling for complete synchronization, while excitatory coupling leads to synchronization with a phase lag. In particular, for inhibitory coupling it was shown that the network syncronizes in multistable clusters of common phase. Wu and Chen (2007, 2009) showed that all-to-all networks with delayed excitatory coupling do not synchronize, either completely or in a weak sense, for sufficiently small delay and coupling strength. In (Wu, Liu, and Chen, 2010) it was shown that the parameter space in systems with excitatory coupling is separated into two regions that support different types of dynamics. The effect of network connectivity to synchronization is numerically studied in (LaMar and Smith, 2010) where it is shown that the proportion of initial conditions that lead to synchronization is an increasing function of the node-degree. Kielblock, Kirst, and Timme (2011) showed that pulses induce the breakdown of order preservation, and demonstrated a system of 2 identical and symmetrically coupled oscillators where the winding numbers of the two oscillators can be different. Klinglmayr and Bettstetter (2012) showed that under self-adjustment assumptions, systems with heterogeneous phases rates and random individual delays would converge to a close-to-synchrony state. Moreover, synchronization has been considered in systems with stochastic features. O’Keeffe, Krapivsky, and Strogatz (2015) studied how small clusters of synchronized oscillators in all-to-all networks coalesce to form larger clusters and obtained exact results for the time-dependent distribution of cluster sizes.

Except for synchronized states more interesting dynamics also manifests in pulse coupled oscillator networks. The existence of unstable attractors has been established, numerically and analytically, in all-to-all pulse coupled oscillator networks with delay, see (Ashwin and Timme, 2005; Broer, Efstathiou, and Subramanian, 2008b; Timme, 2002; Timme, Wolf, and Geisel, 2003). Unstable attractors are fixed points or periodic orbits, which are locally unstable, but have a basin of attraction which is an open subset of the state space. Heteroclinic connections between saddle periodic orbits, such as unstable attractors, have been shown to exist in pulse coupled oscillator networks with delay (Ashwin and Borresen, 2004, 2005; Broer, Efstathiou, and Subramanian, 2008a) and they have been proposed as representations of solutions of computational tasks. Schittler Neves and Timme (2012) showed that complex networks of dynamically connected saddle states are capable of computing arbitrary logic operations by entering into switching sequences in a controlled way. Timme and Wolf (2008) gave an analysis of asymptotic stability for topologically strongly connected PCONs, while Zeitler, Daffertshofer, and Gielen (2009) analyzed the influence of asymmetric coupling and showed that it leads to a smaller bistability range of synchronized states. Zumdieck et al. (2004) numerically showed the existence of long chaotic transients in pulse-coupled oscillator networks. The length of the transients depends on the network connectivity and such transients become prevalent for large networks.

Isochronous dynamics

In this paper we report on a newly observed dynamical behavior of PCONs with delay. Specifically, we show that for appropriate values of the coupling parameters, that is, of the coupling strength $\varepsilon$ and the delay $\tau$, there is a $n\geq 1$-dimensional subset of state space foliated by periodic orbits having the same period. We call the subsets of state space isochronous regions. These periodic orbits are equivalent in a sense we make precise in Definition III.1. Furthermore, the parameter region for which such periodic orbits manifest is an open subset of the parameter space.

This type of observed dynamics in PCON with delay is a special case of isochronous dynamics. One talks of isochronous dynamics when a dynamical system has an open set of initial states that give rise to periodic solutions having the same period. Examples include the one-dimensional harmonic oscillator, any $N$-dimensional harmonic oscillator where the frequencies, $\omega_{1},\dots,\omega_{N}$, satisfy $N-1$ resonance relations, and the restriction of the Kepler problem to any constant energy surface. We refer to (Calogero, 2011) for an extensive review of recent results pertaining to isochronous dynamics in the context of ordinary differential equations and Hamiltonian systems. Nevertheless, such isochronous dynamics have not been previously observed in PCONs, except of course for the trivial case of identical uncoupled oscillators.

A non-trivial example of isochronous dynamics induced by a non-smooth map $g:[0,1]\to[0,1]$ is depicted in Fig. 1. Each point in the middle (red) segment of the graph of $g$, lying along the diagonal, is a fixed point of $g$ and thus such points give isochronous dynamics of period $1$.

Structure of the paper

In Section II we describe the dynamics of PCONs with delay and we review its basic properties. In Section III we first present numerical experiments that show the appearance of $n\geq 1$-dimensional sets of periodic orbits on a surface of section for specific values of the dynamical parameters. Then we define the notion of a isochronous region. In Section IV we discuss in detail one of the isochronous regions in the system. We prove its existence for an open subset of parameter values, describe in detail the dynamics in the region, and determine the stability of the periodic orbits that constitute the region. In Section V we briefly describe other isochronous regions that appear in the system. We conclude the paper in Section VI.

II Dynamics of PCONs with delay

In this section we specify the dynamics of the PCONs with delay that we consider in this paper.

II.1 Mirollo-Strogatz model with delay

We consider a variation of the Mirollo-Strogatz model (Mirollo and Strogatz, 1990; Ernst, Pawelzik, and Geisel, 1995, 1998). The system here is a homogeneous all-to-all network consisting of $N$ pulse coupled oscillators with delayed excitatory interaction. All the oscillators follow the same integrate-and-fire dynamics. Between receiving pulses the state of each oscillator evolves autonomously and its dynamics is smooth. When the $i$-th oscillator reaches the threshold value $x_{i}=1$ its state is reset to $x_{i}=0$. At the same moment the $i$-th oscillator sends a pulse to all other oscillators, $j\neq i$, in the network. The time between the moment an oscillator sends a pulse and the moment the other oscillators receive that pulse is the delay $\tau\geq 0$. When the $i$-th oscillator receives $m$ simultaneous pulses without crossing the threshold value, its state variable jumps to $x_{i}^{\prime}=x_{i}+m\hat{\varepsilon}$. If $x_{i}+m\hat{\varepsilon}\geq 1$, that is, if the oscillator crosses the firing threshold by receiving these pulses, then the new state becomes $x_{i}^{\prime}=0\equiv 1$. The dynamics for each oscillator is thus given by

[TABLE]

Here, $\hat{\varepsilon}=\varepsilon/(N-1)$, where $\varepsilon\geq 0$ is the coupling strength, and $F$ is a positive, decreasing, function $(F>0,\,F^{\prime}<0)$.

To simplify the description of the dynamics we define, following Mirollo and Strogatz (1990), the phases $(\theta_{i})_{i=1}^{N}$ instead of the state variables $(x_{i})_{i=1}^{N}$. The two sets of variables are related through

[TABLE]

where $f:[0,1]\to[0,1]$ is a diffeomorphism fixing the endpoints, that is, $f(0)=0$ and $f(1)=1$. The map $f$ is defined through the requirement that the uncoupled dynamics of each oscillator is given by $\dot{\theta}_{i}=1$. This implies

[TABLE]

and that $f$ is increasing and concave down $(f^{\prime}>0,\,f^{\prime\prime}<0)$. Following Mirollo and Strogatz (1990) we choose

[TABLE]

giving

[TABLE]

Then, in terms of the phases $\theta_{i}$, the dynamics is given by

[TABLE]

The function $H$ is defined by

[TABLE]

and it gives the new phase of an oscillator with phase $\theta$ after it receives a pulse of size $\delta$, ignoring the effect of the threshold.

Typically, one also defines the pulse response function (PRF) $V(\theta,\delta)$ representing the change in phase after receiving a pulse of size $\delta$, ignoring the effect of the threshold. Specifically,

[TABLE]

see Fig. 2b. Note that the function $H$ in Eq. (3) has the property

[TABLE]

implying

[TABLE]

To simplify notation, for fixed value of $\hat{\varepsilon}$, we write

[TABLE]

II.2 Description of the dynamics

In principle, to determine the dynamics of a system with delay $\tau$ for $t\geq 0$ one should know the phases $\theta_{i}(t)$, $i=1,\dots,N$ of the oscillators for all $t\in[-\tau,0]$. This information can be encoded in the phase history function

[TABLE]

In the particular system studied here, this description can be further simplified since it is not all the information about the phases in $[-\tau,0]$ that is necessary to determine the future dynamics. Instead, it is enough to know the phases $\theta_{i}(0)$, $i=1,\dots,N$ at $t=0$ and the firing moments of each oscillator in $[-\tau,0]$, that is, the moments when each oscillator reaches the threshold value.

We denote by $-\sigma_{i}^{(j)}$ the $j$-th firing moment of the $i$-th oscillator in $[-\tau,0]$ and by $\Sigma_{i}=\{\sigma_{i}^{(j)}\}$ the set of all such firings moments. Note that our ordering is

[TABLE]

To simplify notation we also write $\sigma_{i}=\sigma_{i}^{(1)}$ in the case that $\Sigma_{i}$ contains exactly one element. We call the $\sigma_{i}^{(j)}$ firing time distances (FTD) and $\sigma_{i}$ the last firing time distance (LFTD).

Remark II.1**.**

It is shown in Ashwin and Timme (2005) that for sufficiently small values of $\varepsilon$ and $\tau$ the size of the set $\Sigma=\bigcup_{i=1}^{N}\Sigma_{i}$ is bounded for all $t\geq 0$. The parameter region of interest in the present paper is not covered by the explicit estimates given in (Ashwin and Timme, 2005). Nevertheless, for the specific orbits in the isoschronous regions we consider, the size of $\Sigma$ remains bounded for all $t\geq 0$.

The dynamics of the system for $t\geq 0$ can then be determined from the FTD in $[-\tau,0]$ and the phases at $t=0$, i.e., from the set

[TABLE]

When $\Phi$ is a finite set we can ask whether a neighborhood is a finite or infinite dimensional set. Broer, Efstathiou, and Subramanian (2008b) show that, choosing an appropriate metric on the space of phase history functions, a neighborhood of $\Phi$ is finite dimensional. Nevertheless, this local dimension is not constant and is not bounded throughout the state space.

To describe high-dimensional dynamics it is convenient to introduce a Poincaré surface of section. Here we choose the surface $\theta_{N}=0$, see also (Ashwin and Timme, 2005; Broer, Efstathiou, and Subramanian, 2008b). Given a state $\Phi$ with $\theta_{N}=0$ the time evolution of the system produces a new state $\Phi^{\prime}$ when $\theta_{N}$ becomes again [math]. This defines the Poincaré map $\mu:\Phi\to\Phi^{\prime}$. We call the sequence of points $\mu^{j}(\Phi)=\mu\left(\mu^{j-1}(\Phi)\right)$, $j=1,2,\dots$, the Poincaré orbit with initial state $\mu^{0}(\Phi)=\Phi$. We also define a related concept.

Definition II.2** (Phase orbit).**

Consider a Poincaré orbit $\{\mu^{j}(\Phi)\}_{j=0,1,2,\dots}$ and let $\mathrm{pr}_{\theta}$ denote the projection

[TABLE]

Then the phase orbit of $\Phi$ is the sequence $\mathrm{pr}_{\theta}(\mu^{j}(\Phi))$, $j=0,1,2,\dots$.

Remark II.3**.**

Note that the phase orbit gives only a projection of the dynamics to the space of $N-1$ phase variables $(\theta_{1},\dots,\theta_{N-1})$. Since the full dynamics further depends on the firing moments in the time interval $[-\tau,0]$ we cannot define a map $\mathbb{T}^{N-1}\to\mathbb{T}^{N-1}$ that depends only on the phases $(\theta_{1},\dots,\theta_{N-1})$ and fully encodes the dynamics.

A Poincaré orbit $\{\mu^{j}(\Phi)\}_{j=0,1,2,\dots}$ for which $\mu^{j+T_{P}}(\Phi)=\mu_{j}(\Phi)$ for all $j\geq 0$ is called periodic with Poincaré period $T_{P}$. Note that $T_{P}$ is not necessarily the minimal period. By construction, a periodic Poincaré orbit corresponds to a periodic orbit in the full state space for the dynamics with continuous time $t\geq 0$. In particular, let $\Phi(t)$ be the state at time $t\geq 0$ corresponding to a periodic Poincaré orbit. Then there is a time $T$, corresponding to $T_{P}$, such that $\Phi(t+T)=\Phi(t)$ for all $t\geq 0$. We call $T$ the orbit period.

III Isochronous Dynamics

In this paper we consider a pulse coupled oscillator network with $N=3$ oscillators. We show that there is an open region in the paramater space $(\varepsilon,\tau)$ with families of periodic orbits exhibiting intriguing dynamical behavior. In particular, the periodic orbits are not isolated but for each $(\varepsilon,\tau)$ they fill up a $n\geq 2$-dimensional subset in state space, or equivalently, a $n\geq 1$-dimensional subset on the Poincaré surface of section.

III.1 Numerical Experiments

We first report the results of numerical experiments for a pulse coupled $3$-oscillator network with delay with parameters $(\varepsilon,\tau)$. Specifically, we numerically compute the orbits of the system starting from a specific class of initial states $\Phi$ on the Poincaré surface of section $\theta_{3}=0$. These states are defined by scanning the $(\theta_{1},\theta_{2})$-space $\mathbb{T}^{2}$ and setting $\theta_{3}=0$. As we earlier mentioned this information is not sufficient for determining the dynamics of the system and we also need to know the firing time distances. In this computation, for the oscillators $1$ and $2$ we set

[TABLE]

Note that this choice of initial states does not exhaustively cover the phase space due to the restrictions imposed on the FTDs. In particular, we could have also considered initial states with more firing moments in $[-\tau,0]$ but our choice is the simplest natural choice and sufficiently reduces the computational time so as to make the computation feasible while allowing to study the system for different parameter values.

We numerically find that all such orbits are eventually periodic. There is a time $T_{0}$ such that for $t\geq T_{0}$ it holds that $\Phi(t+T)=\Phi(t)$, where $T>0$ is the eventual orbit period. In other words, each initial state converges in finite time to a periodic attractor with period $T$.

In Fig. 3 we show for $(\varepsilon,\tau)=(0.58,0.58)$ the projection of the periodic attractors to the $(\theta_{1},\theta_{2})$-space, that is, we show the phase orbits corresponding to the periodic attractors. The figure shows the existence of periodic orbits with Poincaré periods $T_{P}\in\{3,4,5\}$. Note that we did not find any attractors with Poincaré periods $T_{P}=2$ or $T_{P}\geq 6$ in this computation. Most importantly, we observe that for $(\varepsilon,\tau)=(0.58,0.58)$ the attractors with Poincaré periods $T_{P}\in\{3,4,5\}$ are not isolated. Projections of periodic attractors with $T_{P}=3$ appear to fill one-dimensional sets in the $(\theta_{1},\theta_{2})$-space. Projections of periodic attractors with $T_{P}=4$ or $T_{P}=5$ appear to fill one- and two-dimensional sets. In what follows we analytically study the periodic orbits that we numerically observed. We aim to prove that their projections to the $(\theta_{1},\theta_{2})$-plane fill one- and two-parameter sets and to describe the appearance of these orbits and their properties.

III.2 Definitions

To give a systematic description we classify the periodic orbits into equivalence classes. First, we introduce some notation. Let $O$ be a periodic orbit with period $T>0$ and denote by $P_{i,j}$ the $j$-th pulse received by the $i$-th oscillator in the time interval $[0,T)$. Denote by $n(P_{i,j})$ the multiplicity of the pulse $P_{i,j}$, that is, how many simultaneous pulses correspond to $P_{i,j}$.

Definition III.1**.**

Two periodic orbits $O$ and $O^{\prime}$ are pulse equivalent if they have the same periods $T=T^{\prime}>0$, the sets $\{P_{i,j}\}$ and $\{P^{\prime}_{i,j}\}$ have the same cardinalities, and $n(P_{i,j})=n(P^{\prime}_{i,j})$ for all $i,j$.

We now define an isochronous region. For this we ask that not only the orbit periods in an isochronous region are the same but the stronger condition that the orbits are pulse equivalent.

Definition III.2**.**

A subset $\mathcal{B}$ of the state space is an isochronous region of period $T$ (or Poincaré period $T_{P}$) if

1. (a)

all orbits starting in $\mathcal{B}$ are pulse equivalent with period $T$ (or Poincaré period $T_{P}$), 2. (b)

each orbit starting in $\mathcal{B}$ stays within $\mathcal{B}$, and 3. (c)

there is a homeomorphism $S$ between the space of orbits in $\mathcal{B}$ and an open, connected, subset $\Omega$ of $\mathbb{R}^{k}$, $k\geq 1$.

Remark III.3**.**

$\mathcal{B}$ is required to be invariant under the $\mathbb{R}_{+}$ action induced by the dynamics. This allows to define the space of orbits $\mathcal{B}/\mathbb{R}_{+}$ obtained by reducing $\mathcal{B}$ with respect to the $\mathbb{R}_{+}$ action. Note that the requirement that $\mathcal{B}/\mathbb{R}_{+}$ is connected does not imply that $\mathcal{B}$ is also connected since each periodic orbit in $\mathcal{B}$ may be disconnected. The requirement that $\dim\mathcal{B}/\mathbb{R}_{+}\geq 1$ implies that isolated periodic orbits are excluded.

With these definitions in place, we now turn to the detailed description of one of the isochronous regions that we numerically identified in Section III.

IV The isochronous region IR4

In this section we select one of the numerically observed isochronous regions, describe its periodic orbits, and discuss its existence. In subsection IV.3 we consider the dynamical stability of the periodic orbits. Specifically, we focus on the orbits with $T_{P}=4$ appearing in the lower right corner of Fig. 3b. We denote the corresponding isochronous region by IR4.

IV.1 Description

We have verified, analytically and numerically, that all periodic orbits represented by these points can be parameterized by the firing time distances (FTD) $(\sigma_{1},\sigma_{2},\sigma_{3})$ of the three oscillators. We first prove the following slightly more general result which is also useful for determining the stability of the periodic orbits, see subsection IV.3.

Proposition IV.1** (Dynamics).**

Consider the initial state of the pulse coupled $3$-oscillator network on the Poincaré surface of section $\theta_{3}=0$, determined by the phases $(\theta_{1},\theta_{2},0)$ and the firing time distances $(\{\sigma_{1}\},\{\sigma_{2}\},\{\sigma_{3}\})$ satisfying:

1. (a)

$0<\sigma_{2}<\sigma_{1}<\sigma_{3}<\tau$, 2. (b)

$H_{*}<\theta_{1}+\tau-\sigma_{3}<1$, 3. (c)

$H_{*}<H(\theta_{2}+\tau-\sigma_{3})-\sigma_{1}+\sigma_{3}<1$, 4. (d)

$H_{*}<H(\tau-\sigma_{1})-\sigma_{2}+\sigma_{1}<1$*, * 5. (e)

$H(\sigma_{3}-\sigma_{2})<1$.

Then the dynamics of the system induces the Poincaré map

[TABLE]

Proof.

We use the event sequence representation of the dynamics, see (Broer, Efstathiou, and Subramanian, 2008b). In particular, we denote by $[P,(i_{1},\dots,i_{k}),t]$ a pulse that will be received by the oscillators $i_{1},\dots,i_{k}$ after time $t$. We denote by $[F,i,t]$ the event corresponding to the oscillator $i$ firing after time $t$, further implying that $\theta_{i}=1-t$. The initial condition given by phases $(\theta_{1},\theta_{2},0)$ and firing time distances $(\{\sigma_{1}\},\{\sigma_{2}\},\{\sigma_{3}\})$ corresponds to the event sequence

[TABLE]

Note that we write pulse events separately from fire events, keeping the time ordering in each of the subsets. In particular, this implies that $0<\sigma_{2}<\sigma_{1}<\sigma_{3}<\tau$ and that $0<\theta_{2}<\theta_{1}<1$.

The inequality $\theta_{1}+\tau-\sigma_{3}<1$ implies that the first pulse event will be processed first. Then the next event sequences will be

[TABLE]

Here we used the assumptions that $\theta_{1}+\tau-\sigma_{3}>H_{*}$ and $\theta_{2}+\tau-\sigma_{3}<H_{*}$. The next event sequence is

[TABLE]

The inequality $H(\theta_{2}+\tau-\sigma_{3})-\sigma_{1}+\sigma_{3}<1$ implies again that the first pulse event was processed first and then $H(\theta_{2}+\tau-\sigma_{3})-\sigma_{1}+\sigma_{3}>H_{*}$ that oscillator $2$ fires. Moreover, the assumption $\tau-\sigma_{1}<H_{*}$ ensures that oscillator $3$ does not fire. The next event sequence is

[TABLE]

Here, by the assumptions $H_{*}<H(\tau-\sigma_{1})-\sigma_{2}+\sigma_{1}<1$ and $H(\sigma_{3}-\sigma_{2})+\sigma_{2}<1$, we have

[TABLE]

thus proving the statement. ∎

Let $\Omega_{\varepsilon,\tau}$ be the subset of the $(\sigma_{1},\sigma_{2},\sigma_{3})$-space defined by the relations

[TABLE]

and

[TABLE]

Moreover, define the map

[TABLE]

from $\Omega_{\varepsilon,\tau}$ to the space of initial conditions of the PCON. Then we prove the following statement.

Proposition IV.2**.**

Consider the initial state of the pulse coupled $3$-oscillator network on the Poincaré surface of section $\theta_{3}=0$, given by $S(\bm{\sigma})$ for $\bm{\sigma}\in\Omega_{\varepsilon,\tau}$. Then the map

[TABLE]

has the following properties:

1. (a)

$g(\Omega_{\varepsilon,\tau})=\Omega_{\varepsilon,\tau}$; 2. (b)

$G(S(\bm{\sigma}))=S(g(\bm{\sigma}))$* for all $\bm{\sigma}\in\Omega_{\varepsilon,\tau}$, where $G$ is the Poincaré map (6).*

Proof.

First, one easily checks that if $\bm{\sigma}\in\Omega_{\varepsilon,\tau}$ then $g(\bm{\sigma})\in\Omega_{\varepsilon,\tau}$ and vice versa. Then, note that if $\bm{\sigma}\in\Omega_{\varepsilon,\tau}$ then $S(\bm{\sigma})$ satisfies the conditions of Proposition IV.1. This implies

[TABLE]

∎

Proposition IV.2 shows that $S$ intertwines the map $g$ on $\Omega_{\varepsilon,\tau}$ with the Poincaré map $G$. We then have the following description of the dynamics in $\Omega_{\varepsilon,\tau}$.

Proposition IV.3**.**

The map $g$ on $\Omega_{\varepsilon,\tau}$ has period $4$, that is, $g^{4}(\bm{\sigma})=\bm{\sigma}$ for all $\bm{\sigma}\in\Omega_{\varepsilon,\tau}$. The point $\bm{\sigma}_{*}:=(\tau/2,\tau/4,3\tau/4)\in\Omega_{\varepsilon,\tau}$ is a fixed point of $g$, and points $\bm{\sigma}\in\Omega_{\varepsilon,\tau}$ along the line parameterized by $\bm{\sigma}=\bm{\sigma}_{*}+t(0,1,1)$, $t\in\mathbb{R}$, are period-$2$ points of $g$.

Proof.

The proof of the statement is a straightforward computation. Nevertheless, it is more enlightening to proceed in a different way. Let

[TABLE]

where $\bm{s}=(s_{1},s_{2},s_{3})$. In terms of $\bm{s}$, $g$ becomes the linear map

[TABLE]

where

[TABLE]

Clearly, $\bm{s}=0$ is the only fixed point of $L$. One checks that $L^{2}$ acts as rotation by $\pi$ about the line $\bm{s}=t(0,1,1)$, $t\in\mathbb{R}$. Therefore, $L^{2}$ leaves this line invariant, and $L^{4}$ is the identity. ∎

Remark IV.4**.**

Proposition IV.2 implies that the Poincaré map $G$ has Poincaré period $T_{P}=4$ for each $S(\bm{\sigma})$, $\bm{\sigma}\in\Omega_{\varepsilon,\tau}$. The evolution of the phases of the $3$ oscillators for such orbits is depicted in Fig. 4a and the detailed dynamics is given in Table 1. The set $\Omega_{\varepsilon,\tau}$ also gives rise to periodic orbits with smaller minimal orbit period than $T=3\tau$. In particular, there is a line in $\Omega_{\varepsilon,\tau}$ given by $(\sigma_{1},\sigma_{2},\sigma_{3})=(\tau/2,\sigma_{2},\tau/2+\sigma_{2})$ for which all points give rise to period $T=3\tau/2$ orbits ($T_{P}=2$), see Fig. 4b. One point along this line, having $\sigma_{2}=\tau/4$ gives rise to a period $T=3\tau/4$ orbit ($T_{P}=1$), see Fig. 4c.

IV.2 Existence

Let $\widetilde{\Omega}_{\varepsilon,\tau}=S(\Omega_{\varepsilon,\tau})$ be the embedding of $\Omega_{\varepsilon,\tau}$ in the $(\bm{\theta};\bm{\sigma})$-space. Moreover, let $\mathcal{A}_{\varepsilon,\tau}=\mathrm{pr}_{\theta}(\widetilde{\Omega}_{\varepsilon,\tau})$, where $\mathrm{pr}_{\theta}:\mathbb{R}^{6}\to\mathbb{R}^{2}$ is the projection to the $(\theta_{1},\theta_{2})$-plane. The set $\mathcal{A}_{\varepsilon,\tau}$ is depicted in Fig. 5a for $(\varepsilon,\tau)=(0.58,0.58)$. One can check that $\Omega_{\varepsilon,\tau}$ and $\mathcal{A}_{\varepsilon,\tau}$ have non-empty interior for $(\varepsilon,\tau)=(0.58,0.58)$ and that each point in $\widetilde{\Omega}_{\varepsilon,\tau}$ is the initial condition of a periodic orbit in IR4 of Table 1 with period $T=3\tau$ and Poincaré period $T_{P}=4$ . Therefore, $\widetilde{\Omega}_{\varepsilon,\tau}$ is a periodic plateau.

Proposition IV.5**.**

The isochronous region IR4 of Table 1 exists in the subset of the parameter space $(\varepsilon,\tau)$ given by

[TABLE]

see Fig. 5c.

Proof.

Note that

[TABLE]

This implies that Eq. (10) is a necessary condition for Eq. (7a) to hold. We show that if Eq. (10) holds then $\Omega_{\varepsilon,\tau}$ contains a non-empty open subset. Consider the point

[TABLE]

Then $\bm{\sigma}_{*}\in\Omega_{\varepsilon,\tau}$ if and only if Eq. (10) holds, since in this case we have $F_{k}(\bm{\sigma}_{*};\tau)=H(\tau/2)+\tau/4$, for $k=1,\dots,4$. Therefore, when Eq. (10) holds, $\Omega_{\varepsilon,\tau}\neq\emptyset$. Moreover, when the strict form of Eq. (10) holds, there is an open neighborhood $U$ of $D_{\tau}$ in $\bm{\sigma}$-space such that $U\subset\Omega_{\varepsilon,\tau}$. ∎

Fig. 5d shows the volume of $\Omega_{\varepsilon,\tau}$ for $(\varepsilon,\tau)\in[0,1]^{2}$. The volume is computed using the Mathematica function Volume.

Remark IV.6**.**

If $H(\frac{\tau}{2})+\frac{\tau}{4}=H_{*}$ or $H(\frac{\tau}{2})+\frac{\tau}{4}=1$, the inequalities (7a) are satisfied only by the point $\bm{\sigma}_{*}$. Fig. 6 shows the phases for an orbit starting from the point $\bm{\sigma}_{*}$ when $(\varepsilon,\tau)$ moves outside the region of existence of IR4. In that case the dynamics converges in short time to a stable periodic orbit with $T_{P}=1$.

Fig. 7a compares the analytically obtained $\mathcal{A}_{\varepsilon,\tau}$ for $(\varepsilon,\tau)=(0.58,0.58)$ to the numerical results discussed in subsection III.1. We note that the numerically obtained orbits cover only part of $\mathcal{A}_{\varepsilon,\tau}$. This can be explained by the fact that the space of initial conditions that we scanned in our numerical experiments does not include the periodic orbits in IR4. Some of the orbits in IR4 are periodic attractors for our numerical initial conditions but others are not accessible. This effect is much more pronounced for $(\varepsilon,\tau)=(0.45,0.45)$ as is shown in Fig. 7b. In this case our initial numerical experiments did not reveal the existence of any orbits in IR4. Nevertheless, in this case $\Omega_{\varepsilon,\tau}$ is non-empty and subsequent numerical experiments with different initial conditions allowed us to numerically find orbits in IR4.

IV.3 Stability

In this section we consider the stability of the periodic orbits in IR4. We show that, for $\bm{\sigma}\in\Omega_{\varepsilon,\tau}$, small changes in $\bm{\theta}$ lead to the same periodic orbit while small changes in $\bm{\sigma}$ lead to a nearby periodic orbit in the same pulse equivalence class. In particular, we have the following result.

Proposition IV.7**.**

Let $(\bm{\theta},\bm{\sigma})=S(\bm{\sigma})$, $\bm{\sigma}\in\Omega_{\varepsilon,\tau}$, and denote by $Y(\bm{\sigma})$ the corresponding periodic orbit in IR4. Then for $\Delta\bm{\theta}$ and $\Delta\bm{\sigma}$ sufficiently small, the orbit with initial condition $(\bm{\theta}+\Delta\bm{\theta},\bm{\sigma}+\Delta\bm{\sigma})$ converges in one iteration of the Poincaré map to the periodic orbit $Y(\bm{\sigma}+\Delta\bm{\sigma})$.

Proof.

This statement is a straightforward consequence of Proposition IV.1. Since $\Omega_{\varepsilon,\tau}$ is open in $\mathbb{R}^{3}$, given $\bm{\sigma}\in\Omega_{\varepsilon,\tau}$, there is an open neighborhood $U\ni\bm{\sigma}$ with $U\subseteq\Omega_{\varepsilon,\tau}$. Therefore, for $\Delta\bm{\sigma}$ small enough we have $\bm{\sigma}+\Delta\bm{\sigma}\in\Omega_{\varepsilon,\tau}$. Therefore, $(\bm{\theta}^{\prime},\bm{\sigma}+\Delta\bm{\sigma})=S(\bm{\sigma}+\Delta\bm{\sigma})$ satisfies the conditions of Proposition IV.1. This implies that for sufficiently small $\Delta\bm{\theta}$ we have that $(\bm{\theta}^{\prime}+\Delta\bm{\theta}^{\prime},\bm{\sigma}+\Delta\bm{\sigma})$ also satisfies the conditions of Proposition IV.1 giving convergence to $Y(\bm{\sigma}+\Delta\bm{\sigma})$. Finally, we note that $\Delta\bm{\theta}^{\prime}$ can be made sufficiently small by making $\Delta\bm{\theta}$ sufficiently small because of the continuity of the map $S$. ∎

V Other isochronous regions

The isochronous region IR4 is not the only such region that appears in the system under consideration here. Here we briefly report on two other such regions.

V.1 The isochronous region IR3

The isochronous region IR3 consists of periodic orbits with Poincaré period $T_{P}=3$. For orbits in IR3, two of the oscillators have the same phase. This implies that the projection of orbits in IR3 to the $(\theta_{1},\theta_{2})$-plane lies either on one of the axes or along the diagonal. Let $\Omega_{\varepsilon,\tau}$ be the subset of the $(\sigma_{1},\sigma_{3})$-space defined by the relations

[TABLE]

and

[TABLE]

Then we consider in state space the set $S(\Omega_{\varepsilon,\tau})$ where $S$ is given by

[TABLE]

The region $\Omega_{\varepsilon,\tau}$ and the projection $\mathcal{A}_{\varepsilon,\tau}$ of $S(\Omega_{\varepsilon,\tau})$ on the $(\theta_{1},\theta_{2})$-plane are shown in Fig. 8.

Using similar arguments as in the analysis of IR4 we find that the point

[TABLE]

gives a periodic orbit with Poincaré period $T_{P}=1$. Its phase evolution is shown in Fig. 9. Moreover, we find that this occurs for

[TABLE]

thus giving the subset of the parameter space $(\varepsilon,\tau)$ for which IR3 exists, see Fig. 10a.

V.2 The isochronous region IR5

For the isochronous region IR5, corresponding to periodic orbits with Poincaré period $T_{P}=5$ we consider the subset $\Omega_{\varepsilon,\tau}$ of the $(\sigma_{1},\sigma_{2}^{(1)},\sigma_{2}^{(2)},\sigma_{3})$-space defined by the relations

[TABLE]

Then the set of initial states comprising IR5 is $S(\Omega_{\varepsilon,\tau})$ where $S$ is given by

[TABLE]

The projection of $S(\Omega_{\varepsilon,\tau})$ on the $(\theta_{1},\theta_{2})$-plane is shown in Fig. 12a.

Using similar arguments as in the analysis of IR4 we find that the point

[TABLE]

gives a periodic orbit with Poincaré period $T_{P}=1$. Its phase evolution is shown in Fig. 11. Moreover, we find that this occurs for

[TABLE]

thus giving the subset of the parameter space $(\varepsilon,\tau)$ for which IR5 exists, see Fig. 10b.

VI Conclusions

We have reported the existence of non-trivial isochronous dynamics in pulse coupled oscillator networks with delay. In particular, we have presented numerical evidence for the existence of such isochronous regions and we have proved their existence for a subset of the parameter space $(\varepsilon,\tau)$ with non-empty interior. Moreover, we have described in detail the dynamics and stability of orbits in one of the isochronours regions that we call IR4.

The appearance of isochronous regions in pulse coupled oscillator networks with delays demonstrates the capacity of such systems for generating non-trivial dynamics that one would not, in general, expect for smooth dynamical systems. Of particular interest here is that isochronous dynamics coexists with attracting isolated fixed points and periodic orbits. This may be of interest for applications using heteroclinic connections between saddle periodic orbits as representations of computational tasks (Ashwin and Borresen, 2004, 2005; Schittler Neves and Timme, 2012).

Several questions regarding isochronous regions in pulse coupled oscillator networks with delay remain open. The main questions going forward is whether such dynamics exist for larger numbers of oscillators and whether such dynamics persists in networks with non-identical oscillators or different network structure.

Acknowledgements

This work was completed when P.L., supported by the China Scholarship Council, worked as a visiting PhD student at the University of Groningen. W.L. was supported by the NSFC (Grants no. 11322111 and no. 61273014). K.E. was supported by the NSFC (Grant no. 61502132) and the XJTLU Research Development Fund (no. 12-02-08).