Isochronous Dynamics in Pulse Coupled Oscillator Networks with Delay
Pan Li Wei Lin Konstantinos Efstathiou





P.Li@rug.nl
School of Mathematical Sciences, Centre for Computational Systems Biology of ISTBI, Fudan University, Shanghai 200433, China
Johann Bernoulli Institute for Mathematics and Computer Science,
University of Groningen, P.O. Box 407, 9700 AK, Groningen, The Netherlands


wlin@fudan.edu.cn
School of Mathematical Sciences, Centre for Computational Systems Biology of ISTBI, Fudan University, Shanghai 200433, China
Shanghai Key Laboratory of Contemporary Applied Mathematics, and LMNS, Ministry of Education, China


K.Efstathiou@rug.nl
Johann Bernoulli Institute for Mathematics and Computer Science,
University of Groningen, P.O. Box 407, 9700 AK, Groningen, The Netherlands




  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.




Isochronous Dynamics in Pulse Coupled Oscillator Networks with Delay
    Konstantinos Efstathiou
    December 30, 2023
====================================================================





  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.  




§ 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 <cit.> and the synchronization of fireflies <cit.>. 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 <cit.> who considered excitatory coupling with no delay, <cit.> 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. <cit.> 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 <cit.> 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 <cit.> where it is shown that the proportion of initial conditions that lead to synchronization is an increasing function of the node-degree. <cit.> 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. <cit.> 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. <cit.> 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 <cit.>. 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 <cit.> and they have been proposed as representations of solutions of computational tasks. <cit.> 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. <cit.> gave an analysis of asymptotic stability for topologically strongly connected PCONs, while <cit.> analyzed the influence of asymmetric coupling and showed that it leads to a smaller bistability range of synchronized states. <cit.> 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 ε and the delay τ, there is a n≥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 <ref>. 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, ω_1,…,ω_N, satisfy N-1 resonance relations, and the restriction of the Kepler problem to any constant energy surface. We refer to <cit.> 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] → [0,1] is depicted in <ref>. 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 <ref> we describe the dynamics of PCONs with delay and we review its basic properties. In <ref> we first present numerical experiments that show the appearance of n≥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 <ref> 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 <ref> we briefly describe other isochronous regions that appear in the system. We conclude the paper in <ref>.



§ DYNAMICS OF PCONS WITH DELAY



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



 §.§ Mirollo-Strogatz model with delay


We consider a variation of the Mirollo-Strogatz model <cit.>. 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  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 τ≥ 0. When the i-th oscillator receives m simultaneous pulses without crossing the threshold value, its state variable jumps to x_i' = x_i + m ε̂. If x_i + m ε̂≥ 1, that is, if the oscillator crosses the firing threshold by receiving these pulses, then the new state becomes x_i' = 0 ≡ 1. The dynamics for each oscillator is thus given by

  
    ẋ_i(t)     = F(x_i(t)), 
    
        x_i(t^+)     = 0, if x_i(t)=1, 
    and
        x_i (t)     = min(1, x_i (t^-) + mε̂ ),

  if m other oscillators fired at time t-τ.

Here, ε̂ = ε / (N-1), where ε≥ 0 is the coupling strength, and F is a positive, decreasing, function (F > 0,  F' < 0).

To simplify the description of the dynamics we define, following <cit.>, the phases (θ_i)_i=1^N instead of the state variables (x_i)_i=1^N.
The two sets of variables are related through

    x_i = f(θ_i),

where f:[0,1]→[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 θ̇_i = 1. This implies

    ẋ = F(x) = f'(f^-1(x)) = 1/(f^-1)'(x),

and that f is increasing and concave down (f' > 0,  f” < 0).
Following <cit.> we choose

    F(x) := F_b(x) = e^b - 1/b e^-bx,    b > 0,

giving

    f(θ) := f_b(θ) = 1/bln( 1 + (e^b-1) θ).


Then, in terms of the phases θ_i, the dynamics is given by

  
    θ̇_i(t)     = 1, 
    θ_i(t^+)     = 0, if θ_i(t)=1,
    andθ_i(t)     = min{ 1, H(θ_i(t^-),mε̂) },

  if m other oscillators fired at time t-τ.

The function H is defined by

    H (θ, δ)
      = f^-1(f(θ)+δ)
      = e^b δ θ + e^b δ-1/e^b-1,

and it gives the new phase of an oscillator with phase θ after it receives a pulse of size δ, ignoring the effect of the threshold.

 

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

    V(θ,δ) = H(θ,δ) - θ = (e^bδ -1) θ + e^bδ-1/e^b-1,

see <ref>.
Note that the function H in Eq. (<ref>) has the property

    H(H( θ, δ), δ') = H(θ, δ+δ').

implying

    H(H( θ,  m ε̂),  m' ε̂) = H(θ, (m+m') ε̂).

To simplify notation, for fixed value of ε̂, we write

    H(θ,mε̂) = H_m(θ)   and  H(θ,ε̂) = H_1(θ) = H(θ).




 §.§ Description of the dynamics


In principle, to determine the dynamics of a system with delay τ for t ≥ 0 one should know the phases θ_i(t), i=1,…,N of the oscillators for all t ∈ [-τ,0]. This information can be encoded in the phase history function

    θ: [-τ,0] →𝕋^n : t ↦ (θ_1(t), …, θ_N(t)).


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

We denote by -σ_i^(j) the j-th firing moment of the i-th oscillator in [-τ,0] and by Σ_i = {σ_i^(j)} the set of all such firings moments. Note that our ordering is

    ⋯ < - σ_i^(3) < - σ_i^(2) < - σ_i^(1)≤ 0.

To simplify notation we also write σ_i = σ_i^(1) in the case that Σ_i contains exactly one element. We call the σ_i^(j) firing time distances (FTD) and σ_i the last firing time distance (LFTD).


  It is shown in <cit.> that for sufficiently small values of ε and τ the size of the set Σ = ⋃_i=1^N Σ_i is bounded for all t ≥ 0. The parameter region of interest in the present paper is not covered by the explicit estimates given in <cit.>. Nevertheless, for the specific orbits in the isoschronous regions we consider, the size of Σ remains bounded for all t ≥ 0.


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

    Φ = {{σ_i^(j)}_j, θ_i }_i=1,…,N.

When Φ is a finite set we can ask whether a neighborhood is a finite or infinite dimensional set. <cit.> show that, choosing an appropriate metric on the space of phase history functions, a neighborhood of Φ 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 θ_N = 0, see also <cit.>. Given a state Φ with θ_N = 0 the time evolution of the system produces a new state Φ' when θ_N becomes again 0. This defines the Poincaré map μ: Φ→Φ'. We call the sequence of points μ^j(Φ)=μ( μ^j-1(Φ) ), j=1,2,…, the Poincaré orbit with initial state μ^0(Φ)=Φ. We also define a related concept.


  Consider a Poincaré orbit {μ^j(Φ) }_j=0,1,2,… and let pr_θ denote the projection
  
    Φ = {{σ_i^(j)}_j, θ_i }_i=1,…,N↦{θ_i }_i=1,…,N-1.

  Then the phase orbit of Φ is the sequence pr_θ(μ^j(Φ)), j=0,1,2,….



  Note that the phase orbit gives only a projection of the dynamics to the space of N-1 phase variables (θ_1,…,θ_N-1). Since the full dynamics further depends on the firing moments in the time interval [-τ,0] we cannot define a map 𝕋^N-1→𝕋^N-1 that depends only on the phases (θ_1,…,θ_N-1) and fully encodes the dynamics.


A Poincaré orbit {μ^j(Φ) }_j=0,1,2,… for which μ^j+T_P(Φ) = μ_j(Φ) for all j ≥ 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 ≥ 0. In particular, let Φ(t) be the state at time t ≥ 0 corresponding to a periodic Poincaré orbit. Then there is a time T, corresponding to T_P, such that Φ(t+T) = Φ(t) for all t ≥ 0. We call T the orbit period.




§ 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 (ε,τ) with families of periodic orbits exhibiting intriguing dynamical behavior. In particular, the periodic orbits are not isolated but for each (ε,τ) they fill up a n≥2-dimensional subset in state space, or equivalently, a n≥1-dimensional subset on the Poincaré surface of section.



 §.§ Numerical Experiments



We first report the results of numerical experiments for a pulse coupled 3-oscillator network with delay with parameters (ε, τ). Specifically, we numerically compute the orbits of the system starting from a specific class of initial states Φ on the Poincaré surface of section θ_3 = 0. These states are defined by scanning the (θ_1,θ_2)-space 𝕋^2 and setting θ_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

    Σ_i = {σ_i^(j)}= {θ_i },    if θ_i ≤τ
    ∅,    if θ_i > τ.

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 [-τ,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 ≥ T_0 it holds that Φ(t+T) = Φ(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 <ref> we show for (ε, τ) = (0.58, 0.58) the projection of the periodic attractors to the (θ_1,θ_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 ∈{3,4,5}. Note that we did not find any attractors with Poincaré periods T_P = 2 or T_P ≥ 6 in this computation. Most importantly, we observe that for (ε, τ) = (0.58, 0.58) the attractors with Poincaré periods T_P ∈{3,4,5} are not isolated. Projections of periodic attractors with T_P=3 appear to fill one-dimensional sets in the (θ_1,θ_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 (θ_1,θ_2)-plane fill one- and two-parameter sets and to describe the appearance of these orbits and their properties.



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




  Two periodic orbits O and O' are pulse equivalent if they have the same periods T=T' > 0, the sets { P_i,j} and { P'_i,j} have the same cardinalities, and n(P_i,j) = n(P'_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. 


  A subset ℬ of the state space is an isochronous region of period T (or Poincaré period T_P) if
  
  
  * all orbits starting in ℬ are pulse equivalent with period T (or Poincaré period T_P),
  
  * each orbit starting in ℬ stays within ℬ, and
  
  * there is a homeomorphism S between the space of orbits in ℬ and an open, connected, subset Ω of ℝ^k, k ≥ 1.
  



   ℬ is required to be invariant under the ℝ_+ action induced by the dynamics. This allows to define the space of orbits ℬ / ℝ_+ obtained by reducing ℬ with respect to the ℝ_+ action. Note that the requirement that ℬ / ℝ_+ is connected does not imply that ℬ is also connected since each periodic orbit in ℬ may be disconnected. The requirement that ℬ / ℝ_+ ≥ 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 <ref>.



§ 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 <ref> 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 <ref>. We denote the corresponding isochronous region by IR4.



 §.§ Description


We have verified, analytically and numerically, that all periodic orbits represented by these points can be parameterized by the firing time distances (FTD) (σ_1,σ_2,σ_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 <ref>.


  Consider the initial state of the pulse coupled 3-oscillator network on the Poincaré surface of section θ_3 = 0, determined by the phases (θ_1,θ_2,0) and the firing time distances ({σ_1},{σ_2},{σ_3}) satisfying:
  
  
  * 0 < σ_2 < σ_1 < σ_3 < τ,
  
  * H_* < θ_1 + τ - σ_3 < 1,
  
  * H_* < H(θ_2+τ-σ_3)-σ_1+σ_3 < 1, 
  
  * H_* < H(τ-σ_1)-σ_2+σ_1 < 1,
  
  
  
  * H(σ_3 - σ_2) < 1.
  
  Then the dynamics of the system induces the Poincaré map
  
    G    : (σ_1,σ_2,σ_3;θ_1,θ_2) ↦
          (σ_3-σ_2,σ_1-σ_2,τ-σ_2;H(σ_3-σ_2),σ_1-σ_2).




  We use the event sequence representation of the dynamics, see <cit.>. In particular, we denote by [P,(i_1,…,i_k),t] a pulse that will be received by the oscillators i_1,…,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 θ_i = 1 - t. The initial condition given by phases (θ_1,θ_2,0) and firing time distances ({σ_1},{σ_2},{σ_3}) corresponds to the event sequence
  
    [P,(1,2),τ-σ_3], 
          [P,(2,3),τ-σ_1], 
          [P,(1,3),τ-σ_2],
        
       
          [P,(1,2),τ]; 
          [F,1,1-θ_1], 
          [F,2,1-θ_2], 
          [F,3,1].

  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 < σ_2 < σ_1 < σ_3 < τ and that 0 < θ_2 < θ_1 < 1. 

  The inequality θ_1 + τ - σ_3 < 1 implies that the first pulse event will be processed first. Then the next event sequences will be
  
    1⟶     
          [P,(1,2),0], 
          [P,(2,3),σ_3-σ_1], 
          [P,(1,3),σ_3-σ_2], 
       
          [P,(1,2),σ_3]; 
          [F,1,1-θ_1-τ+σ_3], 
       
          [F,2,1-θ_2-τ+σ_3], 
          [F,3,1-τ+σ_3]
        
    2⟶      [P,(2,3),σ_3-σ_1], 
          [P,(1,3),σ_3-σ_2], 
          [P,(1,2),σ_3]; 
        [F,1,0], 
          [F,2,1-H(θ_2+τ-σ_3)], 
          [F,3,1-τ+σ_3].

  Here we used the assumptions that θ_1 + τ - σ_3 > H_* and θ_2 + τ - σ_3 < H_*. The next event sequence is
  
    3⟶      [P,(2,3),σ_3-σ_1], 
          [P,(1,3),σ_3-σ_2], 
       
          [P,(1,2),σ_3], 
          [P,(2,3),τ]; 
          [F,2,1-H(θ_2+τ-σ_3)], 
       
          [F,3,1-τ+σ_3], 
          [F,1,1]
        
    4⟶      [P,(2,3),0], 
          [P,(1,3),σ_1-σ_2], 
          [P,(1,2),σ_1], 
       
          [P,(2,3),τ+σ_1-σ_3]; 
       
          [F,2,1-H(θ_2+τ-σ_3)+σ_1-σ_3], 
       
          [F,3,1-τ+σ_1], 
          [F,1,1+σ_1-σ_3]
        
    5⟶      [P,(1,3),σ_1-σ_2], 
          [P,(1,2),σ_1], 
       
          [P,(2,3),τ+σ_1-σ_3]; 
          [F,2,0], 
       
          [F,3,1-H(τ-σ_1)], 
          [F,1,1+σ_1-σ_3]
        
    6⟶      [P,(1,3),σ_1-σ_2], 
          [P,(1,2),σ_1], 
       
          [P,(2,3),τ+σ_1-σ_3], 
          [P,(1,3),τ]; 
       
          [F,3,1-H(τ-σ_1)], 
          [F,1,1+σ_1-σ_3], 
          [F,2,1].

  The inequality H(θ_2+τ-σ_3)-σ_1+σ_3 <1 implies again that the first pulse event was processed first and then H(θ_2+τ-σ_3)-σ_1+σ_3 > H_* that oscillator 2 fires. Moreover, the assumption τ - σ_1 < H_* ensures that oscillator 3 does not fire.
  The next event sequence is
  
    7⟶      [P,(1,3),0], 
          [P,(1,2),σ_2], 
          [P,(2,3),τ+σ_2-σ_3], 
       
          [P,(1,3),τ+σ_2-σ_1]; 
       
          [F,3,1-H(τ-σ_1)+σ_2-σ_1], 
       
          [F,1,1+σ_2-σ_3], 
          [F,2,1+σ_2-σ_1]
        
    8⟶      [P,(1,2),σ_2], 
          [P,(2,3),τ+σ_2-σ_3], 
       
          [P,(1,3),τ+σ_2-σ_1]; 
          [F,3,0], 
       
          [F,1,1-H(σ_3-σ_2)], 
          [F,2,1+σ_2-σ_1].

  Here, by the assumptions H_* < H(τ-σ_1)-σ_2+σ_1 < 1 and H(σ_3 - σ_2) + σ_2< 1, we have
  
  
    9⟶      [P,(1,2),σ_2], 
          [P,(2,3),τ+σ_2-σ_3], 
       
          [P,(1,3),τ+σ_2-σ_1], 
          [P,(1,2),τ]
          ; 
        [F,1,1-H(σ_3-σ_2)], 
          [F,2,1+σ_2-σ_1], 
          [F,3,1],

  thus proving the statement.


Let Ω_ε,τ be the subset of the (σ_1,σ_2,σ_3)-space defined by the relations

  
    0 < σ_2 < σ_1 < σ_3 < τ, 
        H_* ≤ F_k(σ; τ) ≤ 1,  k=1,2,3,4,

  where
  
    F_1(σ;τ)    := H(σ_1)+τ-σ_3, 
    
          F_2(σ;τ)    := H(τ-σ_3+σ_2)+σ_3-σ_1, 
    
          F_3(σ;τ)    := H(τ-σ_1)+σ_1-σ_2, 
    
          F_4(σ;τ)    := H(σ_3-σ_2)+σ_2,


and

    H_*    = H_1^-1(1) = e^b-e^bε̂/(e^b-1)e^bε̂.

Moreover, define the map

    S   : (σ_1, σ_2, σ_3)
      ↦  (θ_1, θ_2, θ_3; {σ_1}, {σ_2}, {σ_3}) 
       = (H(σ_1), σ_2, 0; {σ_1}, {σ_2}, {σ_3}),

from Ω_ε,τ to the space of initial conditions of the PCON. Then we prove the following statement.



  Consider the initial state of the pulse coupled 3-oscillator network on the Poincaré surface of section θ_3 = 0, given by S(σ) for σ∈Ω_ε,τ. Then the map
  
    g : ( σ_1, σ_2, σ_3 ) ↦ ( σ_3-σ_2, σ_1-σ_2, τ-σ_2 )

  has the following properties:
  
  
  * g(Ω_ε,τ) = Ω_ε,τ;
  
  * G(S(σ)) = S(g(σ)) for all σ∈Ω_ε,τ, where G is the Poincaré map (<ref>).
  



  First, one easily checks that if σ∈Ω_ε,τ then g(σ) ∈Ω_ε,τ and vice versa. Then, note that if σ∈Ω_ε,τ then S(σ) satisfies the conditions of Proposition <ref>. This implies
  
    G(S(σ))
            = (σ_3-σ_2,σ_1-σ_2,τ-σ_2;  H(σ_3-σ_2),σ_1-σ_2) 
        = S(g(σ)).



Proposition <ref> shows that S intertwines the map g on Ω_ε,τ with the Poincaré map G. We then have the following description of the dynamics in Ω_ε,τ.


  The map g on Ω_ε,τ has period 4, that is, g^4(σ) = σ for all σ∈Ω_ε,τ. The point σ_* := (τ/2,τ/4,3τ/4) ∈Ω_ε,τ is a fixed point of g, and points
  σ∈Ω_ε,τ along the line parameterized by σ = σ_* + t (0,1,1), t ∈ℝ, are period-2 points of g.



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

    σ = σ_* + s,

where s = (s_1,s_2,s_3). In terms of s, g becomes the linear map

    g(s) = L s,

where

    L =
      [  0 -1  1;  1 -1  0;  0 -1  0 ].

Clearly, s=0 is the only fixed point of L. One checks that L^2 acts as rotation by π about the line s = t(0,1,1), t ∈ℝ. Therefore, L^2 leaves this line invariant, and L^4 is the identity.




  Proposition <ref> implies that the Poincaré map G has Poincaré period T_P = 4 for each S(σ), σ∈Ω_ε,τ. The evolution of the phases of the 3 oscillators for such orbits is depicted in <ref> and the detailed dynamics is given in <ref>. The set Ω_ε,τ also gives rise to periodic orbits with smaller minimal orbit period than T = 3τ. In particular, there is a line in Ω_ε,τ given by (σ_1,σ_2,σ_3) = (τ/2,σ_2,τ/2+σ_2) for which all points give rise to period T=3τ/2 orbits (T_P=2), see <ref>. One point along this line, having σ_2 = τ/4 gives rise to a period T=3τ/4 orbit (T_P=1), see <ref>.









 §.§ Existence




Let _ε,τ = S(Ω_ε,τ ) be the embedding of Ω_ε,τ in the (θ;σ)-space. Moreover, let 𝒜_ε,τ = pr_θ(_ε,τ), where pr_θ: ℝ^6 →ℝ^2 is the projection to the (θ_1,θ_2)-plane. The set 𝒜_ε,τ is depicted in <ref> for (ε,τ)=(0.58,0.58). One can check that Ω_ε,τ and 𝒜_ε,τ have non-empty interior for (ε, τ) = (0.58, 0.58) and that each point in _ε,τ is the initial condition of a periodic orbit in IR4 of <ref> with period T = 3τ and Poincaré period T_P = 4 . Therefore, _ε,τ is a periodic plateau.



  The isochronous region IR4 of <ref> exists in the subset of the parameter space (ε, τ) given by
  
    H_* ≤ H( τ/2) + τ/4≤ 1,

  see <ref>.



Note that

    1/4∑_k=1^4 F_k(σ;τ) = H(τ/2) + τ/4.

This implies that Eq. (<ref>) is a necessary condition for Eq. (<ref>) to hold. We show that if Eq. (<ref>) holds then Ω_ε,τ contains a non-empty open subset. Consider the point

    σ_* = (τ/2, τ/4, 3τ/4).

Then σ_* ∈Ω_ε,τ if and only if Eq. (<ref>) holds, since in this case we have F_k(σ_*;τ) = H(τ/2) + τ/4, for k=1,…,4. Therefore, when Eq. (<ref>) holds, Ω_ε,τ∅. Moreover, when the strict form of Eq. (<ref>) holds, there is an open neighborhood U of D_τ in σ-space such that U ⊂Ω_ε,τ.





<ref> shows the volume of Ω_ε,τ for (ε,τ) ∈ [0,1]^2. The volume is computed using the Mathematica function .


  If H(τ/2)+τ/4 = H_* or H(τ/2)+τ/4=1, the inequalities (<ref>) are satisfied only by the point σ_*. <ref> shows the phases for an orbit starting from the point σ_* when (ε,τ) 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.




<ref> compares the analytically obtained 𝒜_ε,τ for (ε,τ) = (0.58,0.58) to the numerical results discussed in <ref>. We note that the numerically obtained orbits cover only part of 𝒜_ε,τ. 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 (ε,τ) = (0.45,0.45) as is shown in <ref>. In this case our initial numerical experiments did not reveal the existence of any orbits in IR4. Nevertheless, in this case Ω_ε,τ is non-empty and subsequent numerical experiments with different initial conditions allowed us to numerically find orbits in IR4.



 §.§ Stability




In this section we consider the stability of the periodic orbits in IR4. We show that, for σ∈Ω_ε,τ, small changes in θ lead to the same periodic orbit while small changes in σ lead to a nearby periodic orbit in the same pulse equivalence class. In particular, we have the following result.


  Let (θ,σ) = S(σ), σ∈Ω_ε,τ, and denote by Y(σ) the corresponding periodic orbit in IR4. Then for Δθ and Δσ sufficiently small, the orbit with initial condition (θ+Δθ,σ+Δσ) converges in one iteration of the Poincaré map to the periodic orbit Y(σ+Δσ).



  This statement is a straightforward consequence of Proposition <ref>. Since Ω_ε,τ is open in ℝ^3, given σ∈Ω_ε,τ, there is an open neighborhood U ∋σ with U ⊆Ω_ε,τ. Therefore, for Δσ small enough we have σ + Δσ∈Ω_ε,τ. Therefore, (θ',σ+Δσ) = S(σ + Δσ) satisfies the conditions of Proposition <ref>. This implies that for sufficiently small Δθ we have that (θ'+Δθ',σ+Δσ) also satisfies the conditions of Proposition <ref> giving convergence to Y(σ+Δσ). Finally, we note that Δθ' can be made sufficiently small by making Δθ sufficiently small because of the continuity of the map S.





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



 §.§ 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 (θ_1,θ_2)-plane lies either on one of the axes or along the diagonal.
Let Ω_ε,τ be the subset of the (σ_1,σ_3)-space defined by the relations

  
    0 < σ_1 < σ_3 < τ, 
        H_* ≤ F_k(σ; τ) ≤ 1,  k=1,2,3, 
        H_**≤ F_k(σ; τ) ≤ 1,  k=4,5,6,

  where
  
    F_1(σ;τ)    := H(σ_2)+τ-σ, 
    
          F_2(σ;τ)    := H(τ-σ_3)+σ_3-σ_2, 
    
          F_3(σ;τ)    := H(σ_3-σ_2)+σ_2, 
    
          F_4(σ;τ)    := σ_3, 
    
          F_5(σ;τ)    := τ-σ_2, 
    
          F_6(σ;τ)    := τ+σ_2-σ_3,


and

    H_**   = H_2^-1(1) = e^b-e^2bε̂/(e^b-1)e^2bε̂.

Then we consider in state space the set S(Ω_ε,τ) where S is given by

    S   : (σ_1, σ_3)
      ↦  (θ_1, θ_2, θ_3; {σ_1}, {σ_2}, {σ_3}) 
       = (σ_1, σ_1, 0; {σ_1}, {σ_1}, {σ_3}).

The region Ω_ε,τ and the projection 𝒜_ε,τ of S(Ω_ε,τ) on the (θ_1,θ_2)-plane are shown in <ref>.







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

    σ_* = ( σ_1, σ_3)
               = (τ/3, 2τ/3),

gives a periodic orbit with Poincaré period T_P=1. Its phase evolution is shown in <ref>. Moreover, we find that this occurs for

    H_* ≤ H( τ/3) + τ/3≤ 1,

thus giving the subset of the parameter space (ε, τ) for which IR3 exists, see <ref>. 




 §.§ The isochronous region IR5







For the isochronous region IR5, corresponding to periodic orbits with Poincaré period T_P=5 we consider the subset Ω_ε,τ of the (σ_1, σ_2^(1), σ_2^(2), σ_3)-space defined by the relations

  
    0 < σ_2^(1) < σ_1 < σ_3 <σ_2^(2) < τ, 
        H_* ≤ F_k(σ; τ) ≤ 1,  k=1,2,3,4,5,

  where
  
    F_1(σ;τ)    := H(σ_2^(1))+τ-σ_3, 
    
          F_2(σ;τ)    := H(τ-σ_2^(2))+σ_2^(2)-σ_1, 
    
          F_3(σ;τ)    := H(σ_2^(2)-σ_3)+σ_3-σ_2^(1), 
    
          F_4(σ;τ)    := H(σ_3-σ_1)+σ_1, 
    
          F_5(σ;τ)    := H(σ_1-σ_2^(1))+σ_2^(1)+τ-σ_2^(2).


Then the set of initial states comprising IR5 is S(Ω_ε,τ) where S is given by

    S     : (σ_1, σ_2^(1),σ_2^(2), σ_3) 
       ↦( θ_1, θ_2, θ_3 ; {σ_1}, {σ_2^(1),σ_2^(2)}, {σ_3}) 
        = ( H(σ_1-σ_2^(1))+σ_2^(1), H(σ_2^(1)), 0 ; 
           {σ_1}, {σ_2^(1),σ_2^(2)}, {σ_3}).

The projection of S(Ω_ε,τ) on the (θ_1,θ_2)-plane is shown in <ref>.

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

    σ_*
      = ( σ_1, σ_2^(1), σ_2^(2), σ_3 )
      = ( 2τ/5, τ/5, 4τ/5, 3τ/5).

gives a periodic orbit with Poincaré period T_P=1. Its phase evolution is shown in <ref>. Moreover, we find that this occurs for

    H_* ≤ H( τ/5) + 2τ/5≤ 1,

thus giving the subset of the parameter space (ε, τ) for which IR5 exists, see <ref>. 





§ 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 (ε,τ) 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 <cit.>.

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




