Synchronization in networks of coupled oscillators using a phase-amplitude framework

Universitat Polit`ecnica de Catalunya Facultat de Matem`atiques i Estad´ıstica Master in Advanced Mathematics and Mathematical Engineering Master’s thesis Synchronization in Networks of Coupled Oscillators using a phase-amplitude framework Laia Pomar Pallar`es Supervised by Gemma Huguet Casades and Marina Vegu´e Llorente January, 2025 I would like to express my gratitude to my advisors, Gemma Huguet and Marina Vegu´e, for giving me the opportunity to work with them, and for their time, support, and encouragement throughout this project. Their guidance and feedback have been essential in shaping the direction of this work. I also want to give special thanks to my friends and family for their support. Abstract Networks of coupled nonlinear oscillators can exhibit a wide range of emergent behaviors depending on the connectivity matrix. In this project, we investigate the dynamics of modular and homogeneous directed networks of oscillators described by the phase-amplitude equations. We employ a recently developed dimension reduction technique (Vegu´e et al., 2023) to analyze various synchronization regimes, such as complete synchronization, cluster synchronization, and desynchronization. Our goal is to adapt this technique to a network of oscillators governed by the phase-amplitude equations and to assess the accuracy of the reduced model through numerical simulations. Keywords Complex networks, dynamical systems, dimensionality reduction, synchronization regimes, phase-amplitude equations, oscillators 1 Contents 1 Introduction 3 2 Statement of the problem 4 2.1 Phase-amplitudereduction .................................. 4 2.1.1 Singleoscillator .................................... 4 2.1.2 Network of coupled oscillators . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 Application to a particular oscillator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3 Dimensionality reduction 11 3.1 Transformation to the complex circle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.2 Spectral reduction in the complex circle . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 4 Results 17 4.1 Analysis of a 2-dimesional network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 4.2 Phase Synchronization Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 4.3 Simulations .......................................... 19 4.3.1 Parameters of the model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4.3.2 Simulations of equilibrium dynamics in the reduced model . . . . . . . . . . . . . 21 4.3.3 Simulations of the complete network . . . . . . . . . . . . . . . . . . . . . . . . . 22 5 Explorations with an alternative model 28 5.1 Alternative coupling functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 5.2 Simulations .......................................... 28 5.2.1 Simulations with homogeneous coupling . . . . . . . . . . . . . . . . . . . . . . . 29 5.2.2 Simulations with the connectivity matrix . . . . . . . . . . . . . . . . . . . . . . . 29 6 Conclusions 33 A Mathematical derivations for the averaging method 35 B Python code used for the simulations 37 2 1. Introduction The study of the dynamics of networks of coupled oscillators is a central topic in various fields of science and engineering, ranging from neuroscience to physics and social dynamics. Synchronization regimes are fundamental to understand the collective behavior of complex systems and have practical applications in fields such as neural synchronization, power grid stability, and biological pattern formation. However, the high-dimensional nature of these networks is a big challenge when performing analytical and computational studies. Constructing a reduced representation of a large complex system remains an open problem, highlighting the value of developing reduction techniques that simplify network dynamics while preserving their essential features. In this project, we focus on applying the spectral reduction introduced in [9] to a network of coupled oscillators represented in phase-amplitude coordinates. The phase-amplitude reduction leverages the inherent properties of limit cycles to describe the dynamics of each oscillator using only two variables: the phase and the amplitude deviation along the least contracting direction. On the other hand, the spectral reduction is a dimensionality-reduction technique for dynamical systems on directed networks organized in clusters, where the nodes within the same cluster have similar activity. We apply these methodologies to a canonical model of an Andronov-Hopf bifurcation, as introduced in [8], to explore whether spectral reduction can accurately approximate the dynamics of networks composed of such oscillators. Our goal is to determine whether the reduced system can faithfully capture the synchronization regimes of the full network. This includes examining the influence of network parameters, such as the strength of general coupling ϵand the phase offset ψ, on the synchronization properties. Additionally, we also begin an investigation into an alternative model with different coupling functions, as presented in [6], setting the stage for future research related to this work. The structure of the project is the following. In Section 2 we introduce the theoretical foundations of oscillators and phase-amplitude reduction, including concepts such as isochrons and isostables, as well as a general form of this parameterization. We also present the equations in Cartesian coordinates of the oscillator presented in [1], and derive its expression in the general form of the phase-amplitude reduction. In Section 3, we apply the dimension-reduction technique that we discussed previously, to the general form of the phase-amplitude parameterization. In Section 4, we study analytical properties of a two-dimensional network of the presented oscillators, introduce a directed homogeneous network divided into two groups, and present some numerical simulations of both the complete network and the reduced two-dimensional system. Finally, in Section 5, we present some open work using an alternative model with the same oscillators but different coupling functions. All the plots shown in this manuscript, without reference, have been developed by the author using Python. The code used for these figures and simulations can be found in the Appendix B. 3 Synchronization in networks of couplet oscillators 2. Statement of the problem 2.1 Phase-amplitude reduction Understanding the dynamics of oscillatory systems often requires modelling complex, high-dimensional systems using simplified approaches. In this section, we introduce phase-amplitude variables, a framework which allows representing the behavior and dynamics of an oscillator x∈Rdusing only two variables that describe the oscillatory phase and deviations from the limit cycle along the slowest contracting direction. The following explanation of phase-amplitude variables is based on [1] and [7]. 2.1.1 Single oscillator We begin by considering a single oscillator that is described by the following system of differential equations ˙x=X(x), (1) where X:Rd→Rdis an analytic vector field. The flow is denoted by ϕt(x), and the system has a T-periodic attracting hyperbolic limit cycle Γ. Our goal is to find a local analytic diffeomorphism K:T× B ⊂ T×Rd−1→Rd (θ,σ)7→ K(θ,σ) = x(2) in such a way that the dynamics of the vector field Xin these new variables is given by ˙ θ=1 T=ω, ˙σ= Λσ, with Λ =    λ1 ... λd−1   , (3) where λ1, ... , λd−1are the characteristic exponents of the periodic orbit. Notice that the variable θrotates at constant speed ωand the variables σicontract at rate λi. It is important to notice that the characteristic exponents λiare negative because we are assuming that the limit cycle Γ is hyperbolic attracting. Moreover, we will suppose that these exponents are real and distinct, that is, λi∈Rand λd−1<... < λ1<0. Using these new variables, we can see that the evolution of the flow ϕt(x) is described by ϕt(K(θ,σ)) = K(θ+ωt,eΛtσ). (4) It is interesting to note that the limit cycle Γ corresponds to the points where σequals zero, i.e., Γ := nx∈Rd|x=K(θ, 0) for θ∈To. (5) Therefore, the function γ(θ) = K(θ, 0) provides a parameterization of the limit cycle Γ in terms of the phase variable θ. As we are considering a hyperbolic limit cycle, for any point pin a neighborhood of the limit cycle Ω contained in its basin of attraction, there exists a point q∈Γ such that lim t→∞|ϕt(p)−ϕt(q)|= 0. (6) 4 We say that pand qhave the same asymptotic phase, and we define the isochron Iθas the set of all points with the same asymptotic phase θ, that is Iθ=np∈Ω|lim t→∞|ϕt(p)−ϕt(γ(θ))|= 0o. (7) The diffeomorphism Ksatisfies the following property ∀x∈ Iθx=K(θ,σ) for some σ∈Rd−1. (8) Moreover, we can define a scalar function Θ that assigns to any point xin Ω, its asymptotic phase θ, that is, Θ: Ω ⊂Rd→T x7→ Θ(x) = θ.(9) With all this, the level curves of Θ correspond to the isochrons Iθof the oscillator Iθ={x∈Ω|Θ(x) = θ}. (10) Analogously, we can also define d−1 scalar functions Σi(i= 1, ..., d−1) that map any point x∈Ω to the amplitude variable σi Σi: Ω ⊂Rd→R x7→ Σi(x) = σi.(11) The level curves of these functions are called isostables Ai σi={x∈Ω|Σi(x) = σi}. (12) To get a better idea of these concepts, in Figure 1, we show an example of some isochrons and isostables of the Andronov-Hopf oscillator presented in [8], which we will study in more detail in the next section, for different values of its parameters. Notice that the vector function (Θ(x), Σ1(x), ... , Σd−1(x)) are the inverse of K, satisfying the following property: K(Θ(x), Σ1(x), ... , Σd−1(x)) = x. (13) So, if we derivative on both sides: Id×d=DK(θ,σ)     ∇Θ(x) ∇Σ1(x) . . . ∇Σd−1(x)      (14) and therefore, DK(θ,σ)−1=     ∇Θ(x) ∇Σ1(x) . . . ∇Σd−1(x)      . (15) 5 Synchronization in networks of couplet oscillators Now for equations (21b) and (21c). ˙σk=λσk+ 2 X m=1 ˜ G(σ) m(zk,zk,σk,hm k), hm k= N X j=1 wkj ˜ Hm(zk,zk,σk,zj,zj,σj). Joining everything, we get the phase-amplitude parameterization of the oscillators on the complex unit circle ˙zk=i2πzkω+ 2 X m=1 i2πzk˜ G(z) m(zk,zk,σk,hm k), (43a) ˙σk=λσk+ 2 X m=1 ˜ G(σ) m(zk,zk,σk,hm k), (43b) hm k= N X j=1 wkj ˜ Hm(zk,zk,σk,zj,zj,σj). (43c) 3.2 Spectral reduction in the complex circle Now that we have defined our system in the complex unit circle (43), we can proceed to apply the spectral reduction. To do this, first we define the following macroscopic observables Zν=X k∈Gν aνkzk, Σν=X k∈Gν aνkσk,Hm ν=X k∈Gν aνkhm kν= 1 ... n, (44) where aν=aνkN k=1 is the reduction vector of the group Gν, and satisfies the following conditions N X k=1 aνk= 1, aνk≥0∀k,aνk= 0 if k/∈Gν. (45) Using equations of our system (43), we derive the exact dynamics of these macroscopic observables ˙ Zν=X k∈Gν aνk˙zk=i2πω X k∈Gν aνkzk+i2π 2 X m=1 X k∈Gν aνkzk˜ G(z) m(zk,zk,σk,hm k) ˙ Zν=i2πωZ+i2π 2 X m=1 X k∈Gν aνkzk˜ G(z) m(zk,zk,σk,hm k), (46) ˙ Σν=X k∈Gν aνk˙σi=λX k∈Gν aνkσk+ 2 X m=1 X k∈Gν aνk˜ G(σ) m(zk,zk,σk,hm k) 12 ˙ Σν=λΣν+ 2 X m=1 X k∈Gν aνk˜ G(σ) m(zk,zk,σk,hm k), (47) Hm ν=X k∈Gν aνkhm k=X k∈Gν aνk N X j=1 wkj ˜ Hm(zk,zk,σk,zj,zj,σj) Hm ν= n X ρ=1 X k∈GνX j∈Gρ aνkwkj ˜ Hm(zk,zk,σk,zj,zj,σj). (48) Our objective is to express this as an ODEs system in terms of the macroscopic variables. However, this cannot be achieved in general. Instead, we approximate the dynamics by using the first-order Taylor series of the functions zk˜ G(z) m(zk,zk,σk,hm k), ˜ G(σ) m(zk,zk,σk,hm k) and ˜ Hm(zk,zk,σk,zj,zj,σj) around the observables in order to get a closed-form system. It is important to note, that this approximation only makes sense when the activities of the nodes within the same group are close to each other, since in that case, the nodes will also be close to their corresponding observables. Let’s start with the functions zk˜ G(z) m(zk,zk,σk,hm k) zk˜ G(z) m(zk,zk,σk,hm k)≈ Zν˜ G(z) m(Zν,Zν, Σν,Hm ν)+(Zν˜ G(z) m,1(Zν,Zν, Σν,Hm ν) + ˜ G(z) m(Zν,Zν, Σν,Hm ν))(zk− Zν) +Zν˜ G(z) m,2(Zν,Zν, Σν,Hm ν)(zk− Zν) + Zν˜ G(z) m,3(Zν,Zν, Σν,Hm ν)(σk−Σν) +Zν˜ G(z) m,4(Zν,Zν, Σν,Hm ν)(hm k− Hm ν), where ˜ G(z) m,l(Zν,Zν, Σν,Hm ν) refers to the derivative of ˜ G(z) m(Zν,Zν, Σν,Hm ν) with respect to the l-th variable. Analogously, ˜ G(σ) m(zk,zk,σk,hm k)≈˜ G(σ) m(Zν,Zν, Σν,Hm ν) + ˜ G(σ) m,1(Zν,Zν, Σν,Hm ν)(zk− Zν) +˜ G(σ) m,2(Zν,Zν, Σν,Hm ν)(zk− Zν) + ˜ G(σ) m,3(Zν,Zν, Σν,Hm ν)(σk−Σν) +˜ G(σ) m,4(Zν,Zν, Σν,Hm ν)(hm k− Hm ν), where ˜ G(σ) m,l(Zν,Zν, Σν,Hm ν) refers to the derivative of ˜ G(σ) m(Zν,Zν, Σν,Hm ν) with respect to the l-th variable. ˜ Hm(zk,zk,σk,zj,zj,σj)≈˜ Hm(Zν,Zν, Σν,Zρ,Zρ, Σρ) + ˜ Hm,1(Zν,Zν, Σν,Zρ,Zρ, Σρ)(zk− Zν) +˜ Hm,2(Zν,Zν, Σν,Zρ,Zρ, Σρ)(zk− Zν) + ˜ Hm,3(Zν,Zν, Σν,Zρ,Zρ, Σρ)(σk−Σν) +˜ Hm,4(Zν,Zν, Σν,Zρ,Zρ, Σρ)(zj− Zρ) + ˜ Hm,5(Zν,Zν, Σν,Zρ,Zρ, Σρ)(zj− Zρ) +˜ Hm,6(Zν,Zν, Σν,Zρ,Zρ, Σρ)(σj−Σρ) where ˜ Hm,l(Zν,Zν, Σν,Zρ,Zρ, Σρ) refers to the derivative of ˜ Hm(Zν,Zν, Σν,Zρ,Zρ, Σρ) with respect to the l-th variable. 13 Synchronization in networks of couplet oscillators Plugging this into equations (46), (47) and (48), respectively, ˙ Zν≈i2πωZν+i2π 2 X m=1 Zν˜ G(z) m(Zν,Zν, Σν,Hm ν), (49) ˙ Σν≈λΣν+ 2 X m=1 ˜ G(σ) m(Zν,Zν, Σν,Hm ν), (50) Hm ν≈ n X ρ=1 "Wνρ ˜ Hm(Zν,Zν, Σν,Zρ,Zρ, Σρ) + ˜ Hm,1(Zν,Zν, Σν,Zρ,Zρ, Σρ) X k∈Gν,j∈Gρ aνkwkj zk− WνρZν  +˜ Hm,2(Zν,Zν, Σν,Zρ,Zρ, Σρ) X k∈Gν,j∈Gρ aνkwkj zk− WνρZν  +˜ Hm,3(Zν,Zν, Σν,Zρ,Zρ, Σρ) X k∈Gν,j∈Gρ aνkwkj σk− WνρΣν  +˜ Hm,4(Zν,Zν, Σν,Zρ,Zρ, Σρ) X k∈Gν,j∈Gρ aνkwkj zj− WνρZρ  +˜ Hm,5(Zν,Zν, Σν,Zρ,Zρ, Σρ) X k∈Gν,j∈Gρ aνkwkj zj− WνρZρ  +˜ Hm,6(Zν,Zν, Σν,Zρ,Zρ, Σρ) X k∈Gν,j∈Gρ aνkwkj σj− WνρΣρ #, (51) where Wνρ := X k∈Gν,j∈Gρ aνkwkj . (52) As we want these equations to depend solely on the macroscopic variables, we impose the following X k∈Gν,j∈Gρ aνkwkj zk=µνρZν,(53a) X k∈Gν,j∈Gρ aνkwkj zk=µνρZν, (53b) X k∈Gν,j∈Gρ aνkwkj σk=µνρΣν, (53c) 14 X k∈Gν,j∈Gρ aνkwkj zj=λνρZρ, (54a) X k∈Gν,j∈Gρ aνkwkj zj=λνρZρ, (54b) X k∈Gν,j∈Gρ aνkwkj σj=λνρΣρ. (54c) These conditions are intuitive, as they show that the weighted sum of terms involving the microscopic activities of kfor all k∈Gνis proportional to the macroscopic observable of this group. Now, instead of (53) and (54), as the microscopic variables are arbitrary, we can simply write, respectively X j∈Gρ wkj aνk=µνρaνk∀k∈Gν, (55) X k∈Gν wkj aνk=λνρaρj∀j∈Gρ. (56) Or, in matrix form Kνρˆ aν=µνρˆ aν, (57) WT νρˆ aν=λνρˆ aρ. (58) Here ˆ aνis a vector which contains the components of the vector aνthat correspond to elements of Gν.Kνρ is the diagonal matrix of size mν×mνwhere (Kνρ)k,k=Pj∈Gρwikj, and Wνρ is a matrix of size mν×mρgiven by Wνρ = (wikjl)k,l, where Gν={i1, ... , imν}and Gρ={j1, ... , jmρ}. From now on, we will refer to equations (57) and (58) as compatibility equations. A comprehensive and formal mathematical explanation of how to solve them approximately can be found in [9]. An interesting observation is that if the functions ˜ Hm(zk,zk,σk,zj,zj,σj) do not depend on zk,zk and σk(or, equivalently, the functions Hm(θk,σk,θj,σj) do not depend on θkand σk), the compatibility equation (57) is not necessary, and we only need to solve (58). This is the case for the oscillator (32) presented in Section 2.2. When 57 and 58 are fulfilled for all ν,ρ∈ {1, ... , n}, we get the approximate dynamics ˙ Zν≈i2πωZν+i2π 2 X m=1 Zν˜ G(z) m(Zν,Zν, Σν,Hm ν), (59a) ˙ Σν≈λΣν+ 2 X m=1 ˜ G(σ) m(Zν,Zν, Σν,Hm ν), (59b) Hm ν≈ n X ρ=1 Wνρ ˜ Hm(Zν,Zν, Σν,Zρ,Zρ, Σρ). (59c) It is important to note that the unknowns in equations (57) and (58) are the reduction vectors ˆaν, as well as the parameters µνρ and λµρ. Therefore, the aim is to determine the best way to weight the 15 Synchronization in networks of couplet oscillators activities within each group, where ”best” refers to ensuring that the dynamics of the observables remain as close as possible to the dynamics of the closed system (59). Additionally, the equations typically do not have an exact solution; the method proposed in the article provides an approximate solution whenever the connection matrices are positive. Finally, we observe that the microscopic dynamics defined in (43) share the same form as the macroscopic dynamics in (59). This indicates that the phases associated with the macroscopic complex observables must follow dynamics analogous to those of the microscopic phases θi. Therefore, we can define macroscopic observables Θνassociated with Zν, which relate to them in the same way as the microscopic variables θirelate to zi. Using these observables, we get the approximate dynamics ˙ Θν≈ω+ 2 X m=1 G(θ) m(Θν, Σν,Hm ν), (60a) ˙ Σν≈λΣν+ 2 X m=1 G(σ) m(Θν, Σν,Hm ν), (60b) Hm ν≈ n X ρ=1 WνρHm(Θν, Σν, Θρ, Σρ). (60c) 16 4. Results We now aim to apply the dimensionality reduction technique to the model introduced in Section 2.2. In this section, we will examine this method in the context of a network comprising Nof these oscillators (32), partitioned into two distinct groups. To achieve this, we first conduct an analytical study of the evolution of a two-dimensional system using averaging techniques. Following this, we perform a series of numerical simulations to compare the synchronization states of the exact and reduced systems, focusing on how these states vary with different parameter values. 4.1 Analysis of a 2-dimesional network Let us consider two weakly coupled oscillators of the form (32). We define φias the phase deviation of the oscillator idue to the interacions from the other oscillators: θi=ωt+φi. (61) By combining (32a) and (61), we can derive an expression for ˙φi ˙φi=ϵ⟨∇Θ(ωt+φi,σi), X j wij F(K(ωt+φj,σj))⟩. (62) For simplicity, abusing of notation, from now on we will consider F(K(ωt+φj,σj)) = F(ωt+φj,σj). Since we are interested in the synchronization of these two oscillators, we define a new variable, χ:= φ2−φ1, which represents the phase difference between them. The dynamics of this variable χwill be the focus of our study, as understanding its behavior is equivalent to analyzing the stability of the various synchronization states between the two oscillators. To do this, we will employ the method of averaging [4], which transforms the system by means of a near-identity change of variable into the form ˙φi=ϵ 2 X j=1 wij T 2 X k=1 ZT 0 ∇Θ(k)(ωt+φi,σi)F(k)(ωt+φj,σj)dt. (63) Defining τ:= ωt+φiand noting that ωt+φj=τ−φi+φj=τ+(φj−φi), we rewrite the expression as ˙φi=ϵ 2 X j=1 wij 2 X k=1 Z1 0 ∇Θ(k)(τ,σi)F(k)(τ+ (φj−φi), σj)dτ. (64) We simplify this to ˙φi=ϵ 2 X j=1 wij H(φj−φi,σi,σj), (65) where H(χ,σi,σj) = 2 X k=1 Z1 0 ∇Θ(k)(τ,σi)F(k)(τ+χ,σj)dτ. (66) Using this, we express the derivative of χas: ˙χ= ˙φ2−˙φ1, (67) 17 Synchronization in networks of couplet oscillators ˙χ=ϵw22H(0, σ2,σ2)−ϵw11H(0, σ1,σ1) + ϵw21H(−χ,σ2,σ1)−ϵw12H(χ,σ1,σ2). (68) As stated earlier, our goal is to find the equilibrium points of this variable. To achieve this, we substitute ∇Θ and Finto H(χ,σi,σj) to obtain a more explicit form. Z1 0 ∇Θ(1)(τ,σi)F(1)(τ+χ,σj)dτ= Z1 0 B(σi)(νsin(ΦK(τ,σi)) + γcos(ΦK(τ,σi)))(K(τ+χ,σj)xcos ψ−K(τ+χ,σj)ysin ψ)dτ, Z1 0 ∇Θ(2)(τ,σi)F(2)(τ+χ,σj)dτ= Z1 0 B(σi)(γsin(ΦK(τ,σi)) −νcos(ΦK(τ,σi)))(K(τ+χ,σj)ycos ψ+K(τ+χ,σj)xsin ψ)dτ. Performing some derivations (provided in Appendix A), we arrive at the following simplified expression H(χ,σi,σj) = 1 2πν Rk(σj) Rk(σi)(νsin D − γcos D) (69) where D= 2πχ +ψ+γln (1 −2ασi) 2ν−γln (1 −2ασj) 2ν, (70) Notice that, when σiequals σj, the function H(χ,σi,σj) becomes independent of these parameters. So we can define the constant: H0:= H(0, σ,σ) = 1 2πν (νsin(ψ)−γcos(ψ)). (71) Therefore, the variable χreaches equilibrium when (w22 −w11)H0+w21H(−χ,σ2,σ1)−w12H(χ,σ1,σ2) = 0. (72) It is important to note that this equilibrium does not depend on the value of ϵ. Using this, we will be able to predict the equilibrium points of our reduced system analytically. This will be explored in detail in Section 4.3.2. 4.2 Phase Synchronization Index When we perform the simulations, we aim to compare the evolution of the exact system with the reduced system. In order to do so, we will use the so-called Synchronization Index (SI), which is a quantitative measure of the synchronization between oscillators. It is defined as [5] r=1 n n X k=1 ei2πθk . (73) . From now on we will consider the case in which n= 2 r=|ei2πθ1+ei2πθ2| 2. (74) 18 θ1 θ2 v 0 3 2π π π 2 Figure 2: Graphical representation of the Synchronization Index (SI). The vectors ei2πθ1and ei2πθ2correspond to the phases of the two oscillators on the unit circle, and their vector sum v=ei2πθ1+ei2πθ2 determines the value of r, calculated as r=1 2|v|. Figure 2 provides a graphical representation of this measure, where the oscillators’ phases θ1and θ2are represented as unit vectors, and the SI is computed using the modulus of their sum. It can be seen that the SI ranges between 0 and 1, where a value of 0 indicates that the two oscillators are in anti-phase, and a value of 1 that they are in phase. It is worth noting that this index depends only on the difference of phase, i.e. θ1−θ2, rather than the two phase values. Let v, be the sum ei2πθ1+ei2πθ2 v=ei2πθ1+ei2πθ2= cos (2πθ1) + isin (2πθ1) + cos (2πθ2) + isin (2πθ2) v= (cos (2πθ1) + cos(2πθ2)) + i(sin (2πθ1) + sin (2πθ2)) |v|2= (cos (2πθ1)+cos (2πθ2))2+(sin (2πθ1)+sin 2πθ2)2= 2+2(cos (2πθ1) cos (2πθ2)+sin (2πθ1) sin (2πθ2)) |v|2= 2(1 + cos (2πθ1−2πθ2)). And using (74), we get the value of the SI in terms of the phase difference r=1 2|v|=r1 + cos (2π(θ1−θ2)) 2. (75) 4.3 Simulations In this section, we present a series of simulations of a network of 100 nodes, comparing the synchronization dynamics observed in the exact model with those in the reduced model. But before diving into the simulations, we first discuss the chosen values for the oscillator parameters and the structure of the connectivity matrix. 19 Synchronization in networks of couplet oscillators 4.3.1 Parameters of the model First, we focus on the parameters chosen for the oscillators. Going back to Section 2.2, we have to set the values of four parameters: α,ν,wand γ. To simplify things, we will choose them in such a way that R= 1 and ω= 1. For this, we set ν=−αand w= 1 −γ. Furthermore, following the paper [1], we will set γ=αa. In this way, the only free parameters are αand a, which we will change along the simulations. These two parameters play an important role in the dynamics of the oscillators. The parameter α controls the strength of the contraction to the limit cycle: for small values of αthe contraction is weak, while for larger values it becomes stronger. On the other hand, the parameter acontrols the transversality of the isochrons to the limit cycle. This can be seen in Figures 3 and 4, respectively. Figure 3 shows that an oscillator with a larger value of αreaches faster the limit cycle (σ= 0), and Figure 4 shows how the curvature of the isochrons increases with the value of the parameter a. Figure 3: Evolution of the amplitude variable σof two oscillators of the form (32) with a= 1 and two different values of α:α= 0.1 (blue), and α= 2 (orange). Figure 4: Isochrons of the studied oscillator (32) for α= 1 and three different values of a:a= 0, a= 2, and a= 5, shown from left to right. Another parameter that has to be set is ψ, which encodes phase information about the connections in the network. In the analysis of the 2-dimensional network that we performed in Section 4.1, we have seen that this parameter affects the points of equilibrium of χ=φ2−φ1. We will study this in more detail in the following section. In addition to the oscillator parameters, we must also specify the details of the network structure. We will work with a network containing 100 nodes, divided evenly into two groups of 50 nodes each. The 20 Figure 5: The left square represents the connectivity matrix, and the right one the in/out-degrees of the nodes in the network. connectivity between nodes is described by a homogeneous connectivity matrix, as well as the parameter ϵ which controls the strength of the general coupling. The matrix is considered homogeneous in the sense that the in-degrees and out-degrees of nodes within the same group exhibit little variability. We can see a representation of this matrix in Figure 5. There are four distinct sections: the top-left section represents the connections within group 1, the top-right section represents the connections from nodes in group 1 to nodes in group 2, the bottom-left section corresponds to the connections from nodes in group 2 to nodes in group 1, and the bottom-right section represents the connections within group 2. Each element wij of the connectivity matrix is equal to 1/Nif there is a connection from node jto node i, or 0 otherwise. However, in practice, for cases where there is no connection, we set wij to a very small value to ensure that the connectivity matrix is strictly positive, as this is required to solve the compatibility equations (57) and (58). For more details on constructing homogeneous networks, see the section Homogeneous Networks in [9]. For the simulations involving the reduced model, where the network is simplified to two dimensions, we utilize the grouped weights as defined in equation (52). 4.3.2 Simulations of equilibrium dynamics in the reduced model In this section, we want to study the equilibrium of χin the reduced model for our parameters. This analysis is performed under the condition σ1=σ2. Since H(χ,σi,σj) does not depend on these parameters when this condition holds, we will write it simply as H(χ). Furthermore, we define G(χ) := w21H(−χ)−w12H(χ). (76) As previously mentioned, the equilibrium points of χdepend on the value of ψ. Figure 6 illustrates the intersection of G(χ) and −H0(w22 −w11) for different values of ψ, representing these equilibrium points. We are also interested in the stability of these equilibrium points, which we study by analyzing the sign of the derivative of G(χ) + H0(w22 −w11) with respect to χ: (G(χ) + H0(w22 −w11))′=G′(χ) = −w21H′(−χ)−w12H′(χ), (77) 21 Synchronization in networks of couplet oscillators 5. Explorations with an alternative model As we discussed in the previous section, the strength of the general coupling ϵdoes not appear to influence the synchronization regimes of the model we studied (26). Motivated by this observation, we aim to explore a model where this parameter plays a significant role. In this subsection, we introduce a different coupling function [6] to the canonical model for the Andronov-Hopf oscillator [8] and outline its parameterization in phase-amplitude variables, which will form the foundation for subsequent simulations and analysis. 5.1 Alternative coupling functions The canonical model for the Andronov-Hopf oscillator [8] with the coupling functions presented in Nicks et al. [6], has the following form, ˙xi=xi−(xi−c2yi)r2+ϵ N X j=1 wij (xi−xj−c1(yj−yi)), (79a) ˙yi=yi−(yi+c2xi)r2+ϵ N X j=1 wij (yi−yj+c1(xj−xi)), (79b) where r=qx2 i+y2 i. We can see that the part of the equations corresponding to the dynamics of a single oscillator is a particular case of the model we studied earlier in (26). Specifically, this occurs when α= 1, ν=−1, w= 0, and γ=−c2, xi−(xi−c2yi)r2=xi(1 −r2)−yi(−c2r2), (80a) yi−(yi+c2xi)r2=yi(1 −r2) + xi(−c2r2). (80b) Thus, we can use the same isomorphism K(θ,σ) defined in (23) to obtain the phase-amplitude parameterization for this model, which will have oscillations of radius R= 1, frequency ω=−c2and Floquet exponent λ=−2. However, the coupling dynamics between the nodes differ from those of the previous model. Consequently, when expressing this system in the general form shown in (21), the only changes pertain to the definitions of the functions Hk, which are now given by [6] H1(θi,σi,θj,σj) = xi−xj−c1(yj−yi), (81) H2(θi,σi,θj,σj) = yi−yj+c1(xj−xi), (82) where (xi,yi) = K(θi,σi) and (xj,yj) = K(θj,σj). 5.2 Simulations Our aim is to find a scenario in which the stability of the synchrony state depends on the general strength of the coupling, ϵ. To this end, we begin by conducting simulations where the connectivity matrix is defined 28 with wij = 1/N, reproducing the scenario studied in Nicks et al. [6], in which they show the stability of the synchrony state for different values of c1,c2and ϵ. After reproducing these results, we extend the analysis by using our network’s connectivity matrix to explore whether there exists a set of parameters c1and c2 for which ϵplays a significant role in regulating the stability of the synchrony state. It is important to note that these explorations remain unfinished and open. The results presented here represent progress towards understanding this phenomenon, but further investigations are needed to fully characterize the interplay between the coupling strength and stability in this alternate model. 5.2.1 Simulations with homogeneous coupling In this section, our objective is to reproduce the study conducted in [6] for a particular set of the parameters c1and c2in order to find a case in which ϵplays a role on the stability of synchrony. Figure 17, presents the findings of this paper about this. The region below the blue curve corresponds to the values of the parameters for which the synchrony state is unstable, and the region above the parameters for which it is stable. Figure 17: Marginal stability curves for synchrony for different parameter values. Extracted from [6]. For our experiments, we set c1=−1.4 and c2= 1.1, and we explored the evolution of the nodes for ϵ= 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, with σand θuniformly distributed in t= 0s. In this case, as the coupling weights were defined as wij =1 Nfor i,j= 1, ... , N, it makes no sense to consider the partition between two groups, and therefore we can’t compute the synchronization index between them. Instead, we present the results by plotting the values of θiand σiat t= 0s, t= 500s and t= 1000s, and also by showing the time evolution of θiin the time window t0= 950s and tf= 1000s. For the simulations in which ϵ≤0.5 the nodes did not reach the synchronization state, as exemplified in Figure 18 for ϵ= 0.5. In contrast, for ϵ= 0.6 and ϵ= 0.7, all nodes ended up synchronized with one another, this is shown in Figure 19 for ϵ= 0.6. This observation is close to, but does not exactly match, the results presented in [6], where the synchronization state for this set of parameters is observed for ϵ≳0.4, as seen in Figure 17. 5.2.2 Simulations with the connectivity matrix Now, we want to find a similar scenario but using the network that we have been studying until now, i.e., setting the weights wij according to the connectivity matrix shown in Figure 5, and considering the partitions into two groups: G1={1, ... , 50}and G2={51, ... , 100}. To explore this, we performed simulations for several sets of parameters c1and c2, varying ϵto investigate the role that the general coupling plays in the synchrony of the system. In Figure 20, we present an 29 Synchronization in networks of couplet oscillators Figure 18: Time evolution of the system (79) with homogeneous weights wij =1 N. Parameters used for the simulation: c1=−1.4, c2= 1.1, and ϵ= 0.5. In t= 0s, the nodes are distributed uniformly. Top: Values of σand θat t= 0s, t= 500s and t= 1000s. Bottom: Evolution of the θvariables in the time window t0= 950s and tf= 1000s. example showing the synchronization index of the exact system for different values of ϵ. We observe that for ϵ≲0.6, the phase difference between the two groups does not appear to reach an equilibrium. The fact that the interval shown on the plot is not just [0, 1] for these cases may just be because the dynamics are too slow and/or the time window being too small. For, ϵ≳0.6 the two groups reach an equilibrium at the synchrony state. In Figure 21, we take a closer look at a case in which the phase difference between the two groups does not stabilize. Similar to what was observed in the other model when ψ= 0, the nodes within the same group present different phase values. This suggests that the reduced system may not be able to capture the dynamics correctly in these cases. However, this part of the work remains open, and to fully characterize the relationship between the coupling strength ϵand the stability of synchrony in this context, as well as to evaluate the performance of the reduced model, a deeper analysis is necessary. This should include additional simulations and also an analytical investigation of the reduced system. 30 Figure 19: Time evolution of the system (79) with homogeneous weights wij =1 N. Parameters used for the simulation: c1=−1.4, c2= 1.1, and ϵ= 0.6. In t= 0s, the nodes are distributed uniformly. Top: Values of σand θat t= 0s, t= 500s and t= 1000s. Bottom: Evolution of the θvariables in the time window t0= 950s and tf= 1000s. Figure 20: Synchronization index for different values of ϵfor the exact dynamics of the system (79) in the time window t0= 950s and tf= 1000s. Parameters used for the simulation: c1=−0.8 and c2= 1.1. 31 Synchronization in networks of couplet oscillators Figure 21: Time evolution of the system for the microscopic variables. The two colors indicate two groups. Parameters used for the simulation: c1=−0.8, c2= 1.1 and ϵ= 0.1. In t= 0, the nodes within the same group have the same value. 32 6. Conclusions In this project, we have studied two main concepts: the phase-amplitude reduction for oscillators, and the application of the spectral reduction presented in Vegu´e et al. [9] to a network of oscillators described by this parameterization. This involved using two macroscopic observables for each group in the network, one representing the phase and the other representing the amplitude. Our objective was to determine whether the reduced system could accurately capture the synchronization regimes of the original network. First, we studied a network of oscillators of the model presented in [8], which corresponds to a canonical model for an Andronov-Hopf bifurcation. Analyzing the results obtained with the numerical simulations, we observed that the reduced system is only able to capture the dynamics of the exact system when the microscopic variables within the same group exhibit similar activities. This is in fact a basic condition in order for the spectral reduction to work effectively, as when the activities within a group are not similar, it does not make sense to approximate all these variables with one macroscopic observable. Additionally, we found that the strength of the general coupling in the network, ϵ, does not affect this synchronization. We have studied this behavior both analytically, for the reduced 2-dimensional system, and through simulations for the complete system. Nevertheless, we have seen that the parameter ψ∈[0, 2π), which encodes phase information about the connections between the nodes, does play a key role in this. Particularly, we have seen that for values of ψapproximately between 4 and 6, the nodes within the same group tend to exhibit similar values of θand σ, and so the reduced system is able to approximate correctly the dynamics of the complete network. Following this, and with the aim of finding a scenario in which the strength of the couplings in the network influences the stability of the synchronization regimes, we started exploring a network with the same oscillators but a different dynamical coupling. We started by replicating part of the study conducted in the paper where the model is presented [6], using homogeneous coupling weights wij . This involved performing simulations for various sets of parameters and comparing our findings with respect to the stability of the synchrony state with the ones shown in the paper. Subsequently, we tried to extend this analysis for our network, employing the connectivity matrix presented in previous simulations. However, this part of the work remains incomplete and requires further exploration. For instance, an analytical study similar to the one conducted in Section 4.1 for a 2-dimensional network with this model would be an interesting path for future research. Additionally, several other future directions could be pursued. One potential extension of the project could be to repeat this study with heterogeneous connectivity matrices, as this type of network is explored in the original paper of the spectral reduction [9]. Moreover, we could study the efficiency of the reduced system on networks with more than two clusters, analyzing the stability of the synchrony and splay states. 33 Synchronization in networks of couplet oscillators References [1] Oriol Castej´on, Antoni Guillamon, and Gemma Huguet. Phase-Amplitude Response Functions for Transient-State Stimuli. J. Math. Neurosci., 3(1):1–26, December 2013. [2] F.C. Hoppensteadt and E.M. Izhikevich. Weakly Connected Neural Networks. Applied Mathematical Sciences. Springer New York, 2012. [3] G. Huguet. Dynamics of coupled oscillators using phase amplitude variables: Efficient algorithms and rigorous results. Personal communication, 2024. [4] Eugene M. Izhikevich. Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. The MIT Press, 2007. [5] Yoshiki Kuramoto. Chemical Oscillations, Waves, and Turbulence. Springer, Berlin, Germany, 1984. [6] R. Nicks, R. Allen, and S. Coombes. Insights into oscillator network dynamics using a phase-isostable framework. Chaos, 34(1), January 2024. [7] Alberto P´erez-Cervera, Tere M-Seara, and Gemma Huguet. Global phase-amplitude description of oscillatory dynamics via the parameterization method. Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(8):083117, 08 2020. [8] Steven H. Strogatz. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry and Engineering. Westview Press, 2000. [9] Marina Vegu´e, Vincent Thibeault, Patrick Desrosiers, and Antoine Allard. Dimension reduction of dynamics on modular and heterogeneous directed networks. PNAS Nexus, 2(5):pgad150, 05 2023. 34 A. Mathematical derivations for the averaging method In Section 4.1, we employed the averaging method to study the synchronization dynamics of a 2-dimensional network. Here, we present in more detail the mathematical derivations performed to obtain Equation (69). Recall that we defined H(χ,σi,σj) = 2 X k=1 Z1 0 ∇Θ(k)(τ,σi)F(k)(τ+χ,σj)dτ. (83) We will start by developing the first term of the summation by substituting the functions ∇Θ(1)(τ,σi) and F(1)(τ+χ,σj) Z1 0 ∇Θ(1)(τ,σi)F(1)(τ+χ,σj)dτ= Z1 0 B(σi)(νsin(ΦK(τ,σi)) + γcos(ΦK(τ,σi)))(K(τ+χ,σj)xcos ψ−K(τ+χ,σj)ysin ψ)dτ= B(σi)Rk(σj)Z1 0 (νsin(ΦK(τ,σi)) + γcos(ΦK(τ,σi)))(cos(ΦK(τ+χ,σj)) cos ψ−sin(ΦK(τ+χ,σj)) sin ψ)dτ. For simplicity, we define A:= ΦK(τ,σi) and B:= ΦK(τ+χ,σj) (This is an abuse of notation, it is important to note that both Aand Bare dependent on τ). Substituting them in the first term we have B(σi)Rk(σj)Z1 0 (νsin A+γcos A)(cos Bcos ψ−sin Bsin ψ)dτ= B(σi)Rk(σj)Z1 0 (νcos ψsin Acos B+γcos ψcos Acos B−νsin ψsin Asin B−γsin ψcos Asin B)dτ= B(σi)Rk(σj) 2Z1 0 cos ψ(ν(sin(A+B) + sin(A−B)) + γ(cos(A+B) + cos(A−B))) −sin ψ(ν(cos(A−B)−cos(A+B)) + γ(sin(A+B)−sin(A−B))) dτ B(σi)Rk(σj) 2Z1 0 sin(A+B)(νcos ψ−γsin ψ) + sin(A−B)(νcos ψ+γsin ψ) + cos(A+B)(γcos ψ+νsin ψ) + cos(A−B)(γcos ψ−νsin ψ)dτ, where A+B= 2π(2τ+χ)−γln(1 −2ασi) + ln(1 −2ασj) 2ν, A−B= 2π(−χ) + γln(1 −2ασi)−ln(1 −2ασj) 2ν. Analogously, for the second term of the summation Z1 0 ∇Θ(2)(τ,σi)F(2)(τ+χ,σj)dτ= Z1 0 B(σi)(γsin(ΦK(τ,σi)) −νcos(ΦK(τ,σi)))(K(τ+χ,σj)ycos ψ+K(τ+χ,σj)xsin ψ)dτ= 35 Synchronization in networks of couplet oscillators B(σi)Rk(σj) 2Z1 0 sin(A+B)(γsin ψ−νcos ψ) + sin(A−B)(γsin ψ+νcos ψ) cos(A+B)(−γcos ψ−νsin ψ) + cos(A−B)(γcos ψ−νsin ψ)dτ. Substituting both terms on equation (83) H(χ,σi,σj) = −Rk(σj) 2πνRk(σi)Z1 0 sin(A−B)(γsin ψ+νcos ψ) + cos(A−B)(γcos ψ−νsin ψ)dτ. Noting that A−Bdoes not depend on τ H(χ,σi,σj) = −Rk(σj) 2πνRk(σi)(sin(A−B)(γsin ψ+νcos ψ) + cos(A−B)(γcos ψ−νsin ψ)) Rk(σj) 2πνRk(σi)(sin(B−A)(γsin ψ+νcos ψ)−cos(B−A)(γcos ψ−νsin ψ)) Rk(σj) 2πνRk(σi)γ 2(−2 cos(B−A+ψ)) + ν 2(2 sin(B−A+ψ)) Rk(σj) 2πνRk(σi)(−γcos(B−A+ψ) + νsin(B−A+ψ)) . Finally, defining D:= B−A+ψ, H(χ,σi,σj) = 1 2πν Rk(σj) Rk(σi)(νsin D − γcos D), (84) where D= 2πχ +ψ+γln (1 −2ασi) 2ν−γln (1 −2ασj) 2ν. (85) 36 B. Python code used for the simulations Here we present the Python code used to perform all the simulations in this project, along with some comments. Imported python libraries 1import numpy as np 2import math 3import cmath 4from scipy.integrate import solve_ivp 5from scipy.optimize import fsolve, root_scalar 6import matplotlib.pyplot as plt Definition of the functions of the model Here we define the functions of the first model we studied as presented in Section 2.2. Functions of the diffeomorphism K(θ,σ). They correspond to equations (23), (24) and (25). 1def diffeo_K(theta, sigma, params): 2# params: parameters of the model (alpha, gamma, nu, psi) 3alpha, gamma, nu, psi = params 4R_k = math.sqrt(alpha/(nu*(2*alpha*sigma - 1))) 5Phi_k = 2*math.pi*theta - (gamma*math.log(1-2*alpha*sigma))/(2*nu) 6return (R_k*math.cos(Phi_k), R_k*math.sin(Phi_k)) 7 8def gradient_theta_1(theta, sigma, params): 9alpha, gamma, nu, psi = params 10 R_k = math.sqrt(alpha/(nu*(2*alpha*sigma - 1))) 11 Phi_k = 2*math.pi*theta - (gamma*math.log(1-2*alpha*sigma))/(2*nu) 12 B = -(2*math.pi*nu*R_k)**(-1) 13 return B*(nu*math.sin(Phi_k)+gamma*math.cos(Phi_k)) 14 15 def gradient_theta_2(theta, sigma, params): 16 alpha, gamma, nu, psi = params 17 R_k = math.sqrt(alpha/(nu*(2*alpha*sigma - 1))) 18 Phi_k = 2*math.pi*theta - (gamma*math.log(1-2*alpha*sigma))/(2*nu) 19 B = -(2*math.pi*nu*R_k)**(-1) 20 return B*(gamma*math.sin(Phi_k)-nu*math.cos(Phi_k)) 21 22 def gradient_sigma_1(theta, sigma, params): 23 alpha, gamma, nu, psi = params 24 R_k = math.sqrt(alpha/(nu*(2*alpha*sigma - 1))) 25 Phi_k = 2*math.pi*theta - (gamma*math.log(1-2*alpha*sigma))/(2*nu) 26 A = (alpha * R_k)/(nu*(2*alpha*sigma-1)) 27 C = -(nu*A)**(-1) 28 return C*math.cos(Phi_k) 29 30 def gradient_sigma_2(theta, sigma, params): 31 alpha, gamma, nu, psi = params 32 R_k = math.sqrt(alpha/(nu*(2*alpha*sigma - 1))) 37 Synchronization in networks of couplet oscillators 8for group_idx, group_slice in enumerate(group_slices): 9sigma_group = (sigma_values[:, group_slice]) 10 size_group = sigma_group.shape[1] 11 for iin range(size_group): 12 plt.plot(timestamps, sigma_group[:, i],color=colors[group_idx]) 13 14 # Formatting the plot 15 plt.title(’Sigmas Time by Groups’) 16 plt.xlabel(’Time’) 17 plt.ylabel(’sigma’) 18 plt.axhline(0, color=’black’, linewidth=0.8, linestyle=’--’)# Line at y=0 for reference 19 plt.ylim(np.min(sigma_values - 0.2),np.max(sigma_values + 0.2)) 20 plt.grid() 21 plt.show() Function to plot the values of θand σat three points in time. These plots show θas the phase variable and σas the radial variable (with an offset so that we can represent negative values of sigma). 1def plot_theta_sigma(theta_values, sigma_values, group_slices, time_plots, timesteps, sigma_interval): 2T,n = theta_values.shape 3theta_scaled = theta_values*2*math.pi 4sigma_min, sigma_max = sigma_interval 5fig, axs = plt.subplots(1, 3, subplot_kw=dict(projection=’polar’), figsize=(15, 5)) 6colors = plt.cm.tab10(range(len(group_slices))) 7offset = 0.2 - sigma_min 8 9sigma_ext = max(-sigma_min,sigma_max) 10 radii_labels = [-sigma_ext, 0, sigma_ext] 11 radii_ticks_all = np.linspace(-sigma_ext, sigma_ext, num=5) + offset 12 13 for i,t_idx in enumerate(time_plots): 14 ax = axs[i] 15 ax.set_title(f"Time = {timesteps[t_idx]}s") 16 17 for group_idx, group_slice in enumerate(group_slices): 18 thetas = theta_scaled[t_idx, group_slice] 19 sigmas = sigma_values[t_idx, group_slice] 20 ax.scatter(thetas, sigmas + offset, label=f’Group {group_idx + 1}’, color=colors[ group_idx]) 21 22 ax.set_ylim(0, np.max(sigma_values+offset)+0.05) # Ensure outer points are not cropped 23 ax.set_yticks(radii_ticks_all) 24 ax.set_yticklabels([]) 25 for rtick in radii_labels: 26 ax.text( 27 0.39, 28 rtick + offset + 0.01, 29 f"{rtick:.2f}", 30 ha=’left’, va=’bottom’, fontsize=10, color="black" 31 ) 32 ax.set_xticklabels([]) 33 ax.spines[’polar’].set_visible(False) 44 34 plt.tight_layout(rect=[0, 0.03, 1, 0.95]) 35 plt.show() This function does the same as the previous one but adding the exact and reduced macroscopic observables. 1def complete_plot_theta_sigma(values, group_slices, time_plots, timesteps, sigma_interval): 2theta_nodes, sigma_nodes, theta_exact_group, sigma_exact_group ,theta_reduced, sigma_reduced = values 3m,n = theta_nodes.shape 4theta_nodes = theta_nodes*2*math.pi 5theta_exact_group = theta_exact_group*2*math.pi 6theta_reduced = theta_reduced*2*math.pi 7sigma_min, sigma_max = sigma_interval 8fig, axs = plt.subplots(1, 3, subplot_kw=dict(projection=’polar’), figsize=(15, 5)) 9colors = plt.cm.tab10(range(len(group_slices))) 10 offset = 0.2 - sigma_min 11 12 sigma_ext = max(-sigma_min,sigma_max) 13 radii_labels = [-sigma_ext, 0, sigma_ext] 14 radii_ticks_all = np.linspace(-sigma_ext, sigma_ext, num=5) + offset 15 16 17 for i,t_idx in enumerate(time_plots): 18 ax = axs[i] 19 ax.set_title(f"Time = {timesteps[t_idx]}s") 20 21 for group_idx, group_slice in enumerate(group_slices): 22 thetas = theta_nodes[t_idx, group_slice] 23 sigmas = sigma_nodes[t_idx, group_slice] 24 ax.scatter(thetas, sigmas + offset, label=f’Group {group_idx + 1}’, color=colors[ group_idx], marker=’+’) 25 26 exact_grouped_theta = theta_exact_group[t_idx][group_idx] 27 exact_grouped_sigma = sigma_exact_group[t_idx][group_idx] 28 reduced_theta = theta_reduced[t_idx][group_idx] 29 reduced_sigma = sigma_reduced[t_idx][group_idx] 30 31 32 ax.scatter( 33 exact_grouped_theta, exact_grouped_sigma + offset, 34 label=f’Exacted grouped (Group {group_idx + 1})’, 35 color=colors[group_idx], marker=’o’, s=100 36 ) 37 ax.scatter( 38 reduced_theta, reduced_sigma + offset, 39 label=f’Reduced (Group {group_idx + 1})’, 40 color=colors[group_idx], marker=’s’, s=100 41 ) 42 43 ax.set_ylim(0, np.max(sigma_nodes+offset)+0.05) # Ensure outer points are not cropped 44 ax.set_yticks(radii_ticks_all) 45 ax.set_yticklabels([]) 46 for rtick in radii_labels: 45 Synchronization in networks of couplet oscillators 47 ax.text( 48 0.39, 49 rtick + offset + 0.01, 50 f"{rtick:.2f}", 51 ha=’left’, va=’bottom’, fontsize=10, color="black" 52 ) 53 ax.set_xticklabels([]) 54 ax.spines[’polar’].set_visible(False) 55 plt.tight_layout(rect=[0, 0.03, 1, 0.95]) 56 plt.show() This function has a similar output as the two previous ones, but it just plots the values of θ 1def plot_theta(theta_values,group_slices, timesteps): 2m,n = theta_values.shape 3theta_scaled = theta_values*2*math.pi 4fig, axs = plt.subplots(1, 3, subplot_kw=dict(projection=’polar’), figsize=(15, 5)) 5colors = plt.cm.tab10(range(len(group_slices))) 6 7 8for iin range(m): 9ax = axs[i] 10 ax.set_title(f"Time = {timesteps[i]}s") 11 12 for group_idx, group_slice in enumerate(group_slices): 13 thetas = theta_scaled[i, group_slice] 14 ax.scatter(thetas, np.ones(len(thetas)), label=f’Group {group_idx + 1}’, color= colors[group_idx]) 15 16 # Format 17 ax.set_yticklabels([]) 18 ax.set_xticklabels([]) 19 ax.spines[’polar’].set_visible(False) 20 plt.tight_layout(rect=[0, 0.03, 1, 0.95]) 21 plt.show() Plotting the orbits Code use to plot the isochrons and isostables of a single oscillator. 1theta_fixed = [0, 0.3, 0.7] 2sigma_min, sigma_max = get_sigma_interval(model_params,5) 3sigma_range = np.linspace(sigma_min, 0.35, 400) 4 5sigma_fixed = [0, -1, 0.3] 6theta_range = np.linspace(0, 1, 400, endpoint=False) 7 8 9colors = plt.cm.tab10(range(6)) 10 plt.figure(figsize=(8, 8)) 11 12 #plot isochrones 46 13 for i, theta in enumerate(theta_fixed): 14 x_vals = [] 15 y_vals = [] 16 for sigma in sigma_range: 17 x, y = diffeo_K(theta, sigma, model_params) 18 x_vals.append(x) 19 y_vals.append(y) 20 plt.plot(x_vals, y_vals, label=f’Isochron for $\\theta$= {theta}’, color=colors[i]) 21 plt.scatter(x_vals, y_vals, color=colors[i], s=10) 22 23 #plot isostables 24 for i, sigma in enumerate(sigma_fixed): 25 x_vals = [] 26 y_vals = [] 27 for theta in theta_range: 28 x, y = diffeo_K(theta, sigma, model_params) 29 x_vals.append(x) 30 y_vals.append(y) 31 plt.plot(x_vals, y_vals, label=f’Isostable for $\\sigma$= {sigma}’, color=colors[i+3], linestyle=’--’) 32 33 plt.xlabel(’x’) 34 plt.ylabel(’y’) 35 plt.title(r’Orbits with Fixed $\theta$’) 36 plt.legend() 37 plt.grid(True) 38 plt.axis(’equal’) 39 plt.show() Averaging analysis Definition of the functions H and G of the averaging analysis. They correspond to equations (69) and (76), respectively. equilibrium target G makes reference to −H0(w22 −w11) as when G(χ0) equals this value, it means that χ0is an equilibrium point. 1def average_H(chi, sigma_i, sigma_j, model_params): 2alpha, gamma, nu, psi = model_params 3aux_i = (gamma*math.log(1-2*alpha*sigma_i))/(2*nu) 4aux_j = (gamma*math.log(1-2*alpha*sigma_j))/(2*nu) 5D = 2*math.pi*chi + psi + aux_i - aux_j 6R_k_i = math.sqrt(alpha/(nu*(2*alpha*sigma_i - 1))) 7R_k_j = math.sqrt(alpha/(nu*(2*alpha*sigma_j - 1))) 8return (1/(2*math.pi*nu))*(R_k_j/R_k_i)*(nu*math.sin(D) - gamma*math.cos(D)) 9 10 def average_G(chi, sigma_i, sigma_j, weights, model_params): 11 return weights[1,0]*average_H(-chi, sigma_j, sigma_i, model_params) - weights[0,1]* average_H(chi, sigma_i, sigma_j, model_params) 12 13 def equilibrium_target_G(sigma_i, sigma_j, weights, model_params): 14 return average_H(0, sigma_j, sigma_j, model_params)*weights[1,1] - average_H(0, sigma_i, sigma_i, model_params)*weights[0,0] 47 Synchronization in networks of couplet oscillators 15 16 def equilibrium_G(chi, sigma_i, sigma_j, weights, model_params): 17 chi = chi[0] if isinstance(chi, (np.ndarray, list)) else chi 18 return average_G(chi, sigma_i, sigma_j, weights, model_params) + equilibrium_target_G( sigma_i, sigma_j, weights, model_params) Definition of the derivatives of H and G, used to assess the stability of the equilibrium points. 1def d_average_H(chi, sigma_i, sigma_j, model_params): 2alpha, gamma, nu, psi = model_params 3aux_i = (gamma*math.log(1-2*alpha*sigma_i))/(2*nu) 4aux_j = (gamma*math.log(1-2*alpha*sigma_j))/(2*nu) 5D = 2*math.pi*chi + psi + aux_i - aux_j 6R_k_i = math.sqrt(alpha/(nu*(2*alpha*sigma_i - 1))) 7R_k_j = math.sqrt(alpha/(nu*(2*alpha*sigma_j - 1))) 8return (1/nu)*(R_k_j/R_k_i)*(nu*math.cos(D) + gamma*math.sin(D)) 9 10 def d_average_G(chi, sigma_i, sigma_j, weights, model_params): 11 return -weights[1,0]*d_average_H(-chi, sigma_j, sigma_i, model_params) - weights[0,1]* d_average_H(chi, sigma_i, sigma_j, model_params) 12 13 14 def stable(root, sigma_i, sigma_j, weights, model_params): 15 return d_average_G(root, sigma_i, sigma_j, weights, model_params) < 0 Code to make the plot in Figure 6 1psi_values = np.linspace(0, 2*math.pi,3, endpoint=False) 2 3chi_values= np.linspace(0, 1, 500) 4G_values_i = [] 5eq_i = [] 6 7for psi_i in psi_values: 8model_params_i = alpha, gamma, nu, psi_i 9G_values = np.array([average_G(chi, sigma_i, sigma_j, grouped_weights, model_params_i) for chi in chi_values]) 10 G_values_i.append(G_values) 11 equilibrium = equilibrium_target_G(sigma_i, sigma_j,grouped_weights, model_params_i) 12 constant_values = np.full(500, equilibrium) 13 eq_i.append(constant_values) 14 15 fig, axes = plt.subplots(1, 3, figsize=(18, 6), sharex=True, sharey=True) 16 17 18 colors = plt.cm.tab10(range(2)) 19 for i, psi in enumerate(psi_values): 20 axes[i].plot(chi_values, G_values_i[i], label=r’G($\chi$)’, color=colors[0]) 21 axes[i].plot(chi_values, eq_i[i], label=r’$-H_0 (w_{22} - w_{11})$’, color=colors[1], linestyle=’--’) 22 axes[i].set_title(f’$\\psi = {psi:.2f}$’) 23 axes[i].set_xlabel(r’$\chi$’) 24 if i == 0: 48 25 axes[i].set_ylabel(’Function Value’) 26 axes[i].grid(True) 27 axes[i].legend() 28 29 30 plt.tight_layout(rect=[0, 0, 1, 0.95]) 31 plt.show() Code to compute the equilibrium points of χfor different ψvalues. 1psi_values = np.arange(0, 2*math.pi,0.1) 2chi_values = np.array([0, 0.5]) 3roots_sigma = [] 4sigma = 0.2 5sigma_i = sigma 6sigma_j = sigma 7min_chi = 1 8significant_psi= -1 9for psi_i in psi_values: 10 roots = [] 11 for chi in chi_values: 12 model_params_i = alpha, gamma, nu, psi_i 13 root = fsolve(equilibrium_G, chi, args=(sigma_i, sigma_j, grouped_weights, model_params_i)) 14 root = root[0]%1 15 if not any(np.isclose(root, existing_root) for existing_root in roots): 16 roots.append(root) 17 if root < min_chi and stable(root, sigma_i, sigma_j, grouped_weights, model_params_i ): 18 min_chi = root 19 significant_psi = psi_i 20 roots_sigma.append(roots) 21 22 print(significant_psi) Code to plot these points based on their stability, as shown in Figure 7 1plt.figure(figsize=(8, 6)) 2 3for i, roots in enumerate(roots_sigma): 4model_params_i = alpha, gamma, nu, psi_values[i] 5stable_roots = [root for root in roots if stable(root, sigma_i, sigma_j, grouped_weights , model_params_i)] 6unstable_roots = [root for root in roots if not stable(root, sigma_i, sigma_j, grouped_weights, model_params_i)] 7 8# Plot stable roots (blue) and unstable roots (red) for each K value 9plt.scatter([psi_values[i]] * len(stable_roots), stable_roots, color=’blue’) 10 plt.scatter([psi_values[i]] * len(unstable_roots), unstable_roots, color=’red’) 11 12 plt.xlabel(’psi’) 13 plt.ylabel(’Roots (y-values)’) 14 plt.title(’Plot of Equilibrium Points Based on Stability’) 49 Synchronization in networks of couplet oscillators 15 16 17 plt.scatter([], [], color=’blue’, label=’Stable Equilibrium’) 18 plt.scatter([], [], color=’red’, label=’Unstable Equilibrium’) 19 20 plt.legend() 21 plt.grid(True) 22 23 plt.show() Synchronization index Functions to compute the synchronization index in a given time window. 1def compute_r(theta1,theta2): 2exp1 = cmath.exp(2*math.pi*1j*theta1) 3exp2 = cmath.exp(2*math.pi*1j*theta2) 4return abs(exp1+exp2)/2 5 6def compute_r_time_window(thetas, time_window): 7r_vector = np.array([compute_r(thetas[time, 0], thetas[time, 1]) for time in time_window]) 8return np.min(r_vector), np.max(r_vector) Code to perform the simulations for different values of ϵ(K in the code). 1R_delta = -1 2sigma_min, sigma_max = get_sigma_interval(model_params,R_delta) 3 4omega = a/(2*math.pi) # frequency 5lambda_ = -2*alpha # contraction value 6 7# micoscopic variables at t = 0 8theta_nodes[0,:] = theta_nodes_0 9sigma_nodes[0,:] = sigma_nodes_0 10 x_nodes_initial = np.concatenate([theta_nodes[0,:], sigma_nodes[0,:]]) 11 12 # macroscopic variables at t = 0 13 (theta_groups[0,:], sigma_groups[0,:]) = get_reduced_observables(reduction_vectors, theta_nodes[0,:], sigma_nodes[0,:]) 14 x_groups_initial = np.concatenate([theta_groups[0,:], sigma_groups[0,:]]) 15 16 K_values = np.arange(0.05,0.95,0.1) 17 18 t_min = T_f-50 19 t_max = T_f 20 time_window = np.where((timesteps >= t_min) & (timesteps <= t_max))[0] 21 22 r_min_exact = np.empty(len(K_values)) 23 r_max_exact = np.empty(len(K_values)) 24 r_min_reduced = np.empty(len(K_values)) 25 r_max_reduced = np.empty(len(K_values)) 26 50 27 for j,Kin enumerate(K_values): 28 print("Computing r for K = ", K) 29 try: 30 # Solve for exact system 31 sol = solve_ivp(derivatives, (0, T_f), x_nodes_initial, t_eval = timesteps, args = ( omega, lambda_, weights*K, model_params,N,R_delta),atol=1e-8, rtol=1e-8) 32 theta_nodes = sol.y[:N, :].T 33 sigma_nodes = sol.y[N:, :].T 34 exact_theta_groups = np.empty((T,2)) 35 exact_sigma_groups = np.empty((T,2)) 36 for iin range(T): 37 (exact_theta_groups[i,:], exact_sigma_groups[i,:]) = get_reduced_observables( reduction_vectors, theta_nodes[i,:], sigma_nodes[i,:]) 38 r_min_exact[j], r_max_exact[j] = compute_r_time_window(exact_theta_groups, time_window) 39 print("r exact: ",r_min_exact[j], r_max_exact[j]) 40 except Exception as e: 41 print(f"Exact system failed for K={K}: {e}") 42 r_min_exact[j] = np.nan 43 r_max_exact[j] = np.nan 44 45 try: 46 # Solve for reduced system 47 sol = solve_ivp(derivatives, (0, T_f), x_groups_initial, t_eval = timesteps, args = (omega, lambda_, grouped_weights*K, model_params,n, R_delta),atol=1e-8, rtol=1e -8) 48 theta_groups = sol.y[:n, :].T 49 sigma_groups = sol.y[n:, :].T 50 r_min_reduced[j], r_max_reduced[j] = compute_r_time_window(theta_groups, time_window ) 51 print("r reduced: ",r_min_reduced[j], r_max_reduced[j]) 52 except Exception as e: 53 print(f"Reduced system failed for K={K}: {e}") 54 r_min_reduced[j] = np.nan 55 r_max_reduced[j] = np.nan Code to plot the results of the previous simulations. 1plt.figure(figsize=(8, 6)) 2 3 4plt.scatter(K_values, r_min_reduced, color=’green’, marker=’o’, label=’r min reduced’) 5plt.scatter(K_values, r_max_reduced, color=’green’, marker=’o’, label=’r max reduced’) 6plt.scatter(K_values, r_min_exact, color=’blue’, marker=’o’, label=’r min exact’) 7plt.scatter(K_values, r_max_exact, color=’blue’, marker=’o’, label=’r max exact’) 8 9for iin range(len(K_values)): 10 plt.plot([K_values[i], K_values[i]], [r_min_exact[i], r_max_exact[i]], color=’blue’, lw =1, linestyle=’--’) 11 plt.plot([K_values[i], K_values[i]], [r_min_reduced[i], r_max_reduced[i]], color=’green’ , lw=1, linestyle=’--’) 12 13 plt.xlabel(r’$\epsilon$’) 51 Synchronization in networks of couplet oscillators 14 plt.ylabel(’r’) 15 plt.title(’’) 16 plt.legend() 17 plt.grid(True) 18 plt.ylim(-0.1, 1.1) 19 20 plt.show() 52