Phase Response Curves of Coupled Oscillators
aa r X i v : . [ q - b i o . N C ] S e p Phase Response Curves of Coupled Oscillators
Tae-Wook Ko ∗ and G. Bard Ermentrout † Department of Mathematics, University of Pittsburgh, Pittsburgh, Pennsylvania 15260, USA (Dated: October 30, 2018)Many real oscillators are coupled to other oscillators and the coupling can affect the response ofthe oscillators to stimuli. We investigate phase response curves (PRCs) of coupled oscillators. ThePRCs for two weakly coupled phase-locked oscillators are analytically obtained in terms of the PRCfor uncoupled oscillators and the coupling function of the system. Through simulation and analyticmethods, the PRCs for globally coupled oscillators are also discussed.
PACS numbers: 05.45.Xt, 89.75.-k, 87.19.La
Many systems in physics, chemistry and biology aremodeled as interacting nonlinear oscillators [1, 2, 3, 4, 5,6]. One of the easiest ways to characterize an oscillator isits phase response curve (PRC)[3, 4, 5, 6, 7]. The PRC isdefined as the steady phase shift of an oscillation relativeto the unperturbed oscillation as a function of the timingof perturbation to the oscillator. It provides a usefulinformation for understanding the oscillator’s behaviorwhen the oscillator is subjected to external stimuli orsignals from other oscillators.In most of previous studies, the PRC is obtained whenthe oscillator is isolated from other oscillators [3, 4, 5, 7].However, many oscillators in real systems are coupled toothers when they are under the influence of external stim-uli, and the coupling can affect the response of the oscil-lators. To better understand the dynamics of oscillatorssuch as the response of neuronal population to signalsfrom other brain region [5] or to controlling stimulations[6], it is necessary to study how the coupling changesthe PRCs. This study can also give insights into thephase response of a giant oscillator (for example, circa-dian rhythm generators [3]) composed of many individualoscillators [8]. In this letter, we study the PRC of cou-pled oscillators using the average phase of the system andthe relative phases between the oscillators comprising thesystem. The PRC is shown to depend on the PRC of theisolated oscillator, the nature of the coupling, and therelative phases between the oscillators. For some cases,the PRCs are analytically obtained. Our approach differsfrom that of Ref. [8] in that we analytically approximatethe PRC while they require the numerical evaluation ofthe adjoint of a certain linear operator.If coupling between a network of oscillators is suffi-ciently “weak”, the possibly high-dimensional system canbe reduced to a network of coupled phase models [2, 4, 5].In the following we exploit this fact and restrict our anal-ysis to coupled phase models. Consider, first, two weaklycoupled phase-locked oscillators subjected to a common ∗ Electronic address: [email protected] † Electronic address: [email protected] perturbation characterized by their individual PRC:˙ θ = ω + KH ( θ − θ ) + Z ( θ ) · Aδ ( t − t ) , (1)˙ θ = ω + KH ( θ − θ ) + Z ( θ ) · Aδ ( t − t ) , (2)where θ i ( t ) is the phase of oscillator i at time t , ω i is thenatural frequency of the oscillator i and K ( ≥
0) is thecoupling strength. H ( θ ) is the coupling function obtainedby the phase reduction [2, 4, 5]. Aδ ( t − t ) denotes aDirac delta impulse with amplitude A at time t which issufficiently large so that the perturbing impulse is appliedafter the system reaches a steady state. Z ( θ ) is the PRCfor uncoupled oscillator obtained using an impulse withunit amplitude. Without coupling ( K = 0), the impulsecauses steady phase shift AZ ( θ i ( t )) for oscillator i .In the presence of coupling ( K = 0), if the oscillatorsare locked with nonzero phase difference, or the inputamplitudes are different, then the input impulse gener-ally causes nonidentical phase changes to the oscillators.Thus, the system transiently deviates from the lockedstate and then returns to the state. The coupling can af-fect the phase shift which the oscillation of the recoveredstate can have relative to the unperturbed oscillation. Wewish to determine the PRC of the coupled oscillators, inother words, how the phase shift depends on the phaseat t of the perturbation.To analyze the dynamics, we convert Eqs. (1) and(2) into those for the average phase Φ ≡ θ + θ and therelative phase φ ≡ θ − θ .˙Φ = ¯ ω + KH e ( φ ) + Z av ( θ , θ ) · Aδ ( t − t ) , (3)˙ φ = ∆ ω − KH o ( φ ) + Z d ( θ , θ ) · Aδ ( t − t ) , (4)where ¯ ω = ω + ω , ∆ ω = ω − ω , Z av ( θ , θ ) = Z ( θ )+ Z ( θ )2 , Z d ( θ , θ ) = Z ( θ ) − Z ( θ ), H e ( φ ) = H ( φ )+ H ( − φ )2 , and H o ( φ ) = H ( φ ) − H ( − φ )2 .For simplicity, let us assume that the system has onestable locked state with φ = φ satisfying 0 = ∆ ω − KH o ( φ ) and H o ′ ( φ ) >
0. The phase of each oscillatorcan be written as θ = Φ + φ/ θ = Φ − φ/
2. Let usdenote the phase shift in a phase, for example θ , relativeto the unperturbed oscillation by ∆ θ . We can see thatthe phase shift for the oscillator 1 is given by∆ θ = ∆Φ + ∆ φ/ . (5) -1-0.8-0.6-0.4-0.2 0 0.2 0.4 0.6 0.8 1 -0.4-0.2 0 0.2 0.4 0.6 0.8 1 -1-0.8-0.6-0.4-0.2 0 0.2 0.4 0.6 0.8 1 ZAZ c1 simtheory -1-0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 θ θ θ θ (a)(c) (b)(d) PR C PR C PR C PR C FIG. 1: (Color online) PRCs with odd coupling functions. ω = π/ ω and ω = π/
2. ∆ ω = 0 . A = 1.(a) K = 1: φ ≈ . K = 0 . φ ≈ . K = 0 . φ ≈ . H ( θ ) = sin θ and Z ( θ ) = − sin θ . (d) H ( θ ) = sin θ − . θ ), Z ( θ ) = − [sin( θ + 0 . π ) − sin(0 . π )] / [1 + sin(0 . π )], and K = 0 . φ ≈ . When H ( θ ) is an odd function, the average phase Φevolves with a constant frequency ¯ ω before and after theimpulse (Eq. (3)). Thus, ∆Φ = AZ av ( θ ( t ) , θ ( t )).When the relative phase remains in the basin of attrac-tion of the original relative phase right after the impulse, φ approaches the original relative phase. Otherwise, therelative phase moves to another stable value (called walk-through). Thus, ∆ φ = φ f − φ , where φ f is the stablevalue of φ reached after the impulse. Note that even φ = φ and φ = φ + 2 π give different results. Therefore,the PRC of the oscillator 1 in the coupled cases is givenby Z c ( θ ) = AZ av ( θ , θ − φ ) + ( φ f − φ ) / . (6)We simulate Eqs. (1) and (2) using Euler method withtime step ∆ t = 0 .
01. We measure the steady phase shiftdue to the impulse relative to the unperturbed activity.The PRC is given by this phase shift as a function of thephase at which the impulse is applied.Figure 1 shows Z c ( θ ) with odd coupling functions.The prediction from the theory (black solid curves)matches very well with the simulation results (symbols).With larger values of ∆ ω and/or smaller values of K , theoscillators are locked with larger φ . In Figs. 1(a), (b),and (c) with H ( θ ) = sin θ and Z ( θ ) = − sin θ , we showthe PRC for different values of the coupling strength K .When φ is very small, the PRC of coupled oscillatorsis very close to that of uncoupled oscillators as expected(Fig. 1(a)). In this case, φ goes to the original value φ after the impulse. In Fig. 1(b), with the larger φ ,the PRC of the coupled oscillator becomes significantlydifferent from that of uncoupled oscillators. When theimpulse can kick the system out of the basin of the sta-ble locked state with φ , the system goes through phasewalk through. If the system has a stable fixed point with -0.4-0.2 0 0.2 0.4 0.6 0.8 1 1.2 ZAZAZ c1av simtheory -0.15-0.1-0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 -0.4-0.2 0 0.2 0.4 0.6 0.8 1 -0.1 0 0.2 0.4 -0.1 0 0.2 0.4 θ θ θ θ (a)(c) (b)(d) PR C PR C PR C PR C FIG. 2: (Color online) PRCs with non-odd couplingfunctions. For (a)-(c), ω = π/ ω , ω = π/ ω = 0 .
4. (a) H ( θ ) = sin( θ + 0 . π ), Z ( θ ) = − [sin( θ + 0 . π ) − sin(0 . π )] / [1 + sin(0 . π )], and K = 1: φ ≈ . A = 1. (b) H ( θ ) = sin( θ − . π ) + 0 . θ − . π ), Z ( θ ) = − [sin( θ + 0 . π ) − sin(0 . π )] / [1 + sin(0 . π )],and K = 0 . φ ≈ . A = 0 .
2. (c) H ( θ ), Z ( θ ), and K are the same as in (b) : φ ≈ . A = 0 .
5. (d) gap-junction coupled Morris-Lecar oscillators [4]: I ext, = I + ∆ I and I ext, = I − ∆ I with ∆ I = 0 . A V = 40. (top) type Icase, I = 50 and g syn = 0 . φ ≈ . I = 94 and g syn = 0 . φ ≈ . φ and an unstable fixed point φ u in [0 , π ) as in thecase with H ( θ ) = sin θ for ∆ ω < K , φ u has the roleof basin boundary and φ goes to φ f = φ + 2 π when φ + AZ ( θ ( t )) − AZ ( θ ( t )) > φ u . This type of changesin φ causes the discontinuity shown in the PRC of Fig.1(c). We show similar results for a coupling function withhigher order Fourier terms and an asymmetric PRC (Fig.1(d)).When H ( θ ) is not an odd function, the even part of H affects the dynamics of Φ and thus the phase shift∆Φ through Eq. (3). Finding the PRC in the ana-lytic form is not possible for these cases since we haveto solve equation (4) for general initial data. Instead,we can get an approximation of the PRC in the limitof small changes in φ . Let φ = φ + q with | q | ≪ q ≈ q e − KH o ′ ( φ )( t − t ) for t > t where q is the changein φ right after the impulse: q = AZ d ( θ ( t ) , θ ( t )).As t → ∞ , φ returns to φ . Thus, ∆ φ = 0. Thephase shift ∆Φ is given by ∆Φ = AZ av ( θ ( t ) , θ ( t )) + R ∞ t K [ H e ( φ ) − H e ( φ )] dt ≈ AZ av ( θ ( t ) , θ ( t )) + H e ′ ( φ )2 H o ′ ( φ ) q , where we use H e ( φ ) − H e ( φ ) ≈ H e ′ ( φ ) q ( t ).Therefore, the PRC of the oscillator 1 is Z c ( θ ) ≈ AZ av ( θ , θ − φ )+ H e ′ ( φ )2 H o ′ ( φ ) AZ d ( θ , θ − φ ) . (7)Figure 2 shows Z c ( θ ) with non-odd coupling func-tions. In Fig. 2(a), we show the PRC with the sim-ple type of H . While AZ and AZ av are similar, theobtained PRC for the coupled oscillator is significantlydifferent from them. The curve from Eq. (7) fits wellwith simulation results for the entire range of θ . Fig-ure 2(b) shows the results with a H function with higherorder terms. We use small A = 0 . A = 0 .
5. With the larger A , the theorymismatches significantly for a range of phases, but stillgives a relatively similar shape to the simulations. Theoverall matching is due to the fact that q = 0 at somephases satisfying Z ( θ ) = Z ( θ − φ ) and around thosephases the theory fits well with simulation results. Fig-ure 2(d) shows the PRC of gap-junction coupled Morris-Lecar oscillators with slightly different injection currents[4]: C ˙ V i = − I ( V i , w i )+ I ext,i + g syn ( V j − V i )+ A V δ ( t − t )with j = 2 , i = 1 ,
2. The details are in Ref.[4]. The system is simulated using the 4th-order Runge-Kutta method. For type I and type II cases, the theorygives good fitting with the simulation results with weakstimulus.Next, we want to understand PRCs for oscillators cou-pled to many other oscillators. We study the case withglobally coupled oscillators: For i = 1 , , ..., N, ˙ θ i = ω i + KN N X j =1 H ( θ j − θ i ) + Z ( θ i ) · Aδ ( t − t ) , (8)where N is the total number of oscillators and others areas defined in the two oscillator system.We introduce similar variables as in two coupled oscil-lators: the average of the phases Φ ≡ N P Nj =1 θ j and thephase φ i ≡ θ i − θ M of oscillator i relative to the phaseof oscillator M , where the subscript M denotes the os-cillator which has the average frequency ¯ ω . From thedefinitions of Φ and φ , we obtain θ M = Φ − N P Nj =1 φ j .Because the PRC for other oscillators can be treated sim-ilarly and oscillator M follows closely to the collectivebehavior of the system, we focus on the PRC Z c ( θ M ) ofoscillator M . As in the case of two coupled oscillators,we get Z c ( θ M ) = ∆Φ − h ∆ φ i with h ∆ φ i ≡ N N X j =1 ∆ φ j . (9)The equation for Φ is˙Φ = ¯ ω + KN N X i,j =1 H e ( φ j − φ i ) + Z av · Aδ ( t − t ) , (10)where Z av ( θ , . . . , θ N ) ≡ N P Ni =1 Z ( θ i ( t )).For an odd function H ( θ ), ∆Φ = AZ av . But for anon-odd function H ( θ ), the second term contributes to∆Φ and it is not easy to calculate ∆Φ analytically.Let us consider fully locked states first. For a fullylocked state with relative phases φ i , the system returnsto the locked state after the stimulation and the relative phase φ i can be changed to the equivalent phase φ i +2 n i π , where n i is an integer. Thus, ∆ φ i = 2 n i π .With H ( θ ) = sin( θ + β ), which is a good approximationfor many general coupling functions, Eq. (8) becomes˙ θ i = ω i + KR sin(Θ − θ i + β ) + Z ( θ i ) · Aδ ( t − t ) , (11)where R and Θ are the order parameter and the corre-sponding collective phase respectively defined by Re i Θ ≡ N P Nj =1 e iθ j . In the frame rotating with the synchroniza-tion frequency Ω, the equation becomes˙ ψ i = ω i − Ω + KR sin( ˜Θ − ψ i + β )+ Z ( θ i ) · Aδ ( t − t ) , (12)where ψ i ≡ θ i − (Ω t + Θ ), ˜Θ ≡ Θ − (Ω t + Θ ) and theconstant Θ is chosen such that the stationary value of ˜Θbefore the impulse is zero. Let R denote the stationaryvalue of R . We can analyze the stationary state of thesystem, using self-consistency argument and find R andΩ [2, 9].For fully locked states or partially locked states wherethe oscillator M is locked with a locking phase ψ M ∗ andthe oscillators form a stationary distribution relative tothe frame rotating with Ω, Z av = N P Ni =1 Z ( θ M + ψ i − ψ M ∗ ), where ψ M ∗ = sin − ( ¯ ω − Ω KR ) + β . For Z ( θ ) = a + a sin( θ + ξ ), using R = N P Nj =1 e i ψ j , we obtain Z av = a + a R sin( θ M + ξ − ψ M ∗ ) . (13)Note that since R <
1, the magnitude of the sinusoidalpart of Z av is smaller than that of Z unless synchrony isperfect.Figures 3(a)-(e) show results with H ( θ ) = sin θ and auniform distribution for the frequencies of the oscillators.We use Z ( θ ) = − sin θ for (a)-(d), and an asymmetric Z ( θ ) for (e). With the given coupling strength, the sys-tem shows a fully locked state (Fig. 3(a)). Figures 3(b)and (c) show the PRCs for different values of A . Withweak stimulation (Fig. 3(b)), all φ i return to the unper-turbed values (∆ φ i = 0 for all i , Fig. 3(d)) and the PRCis shown to be contributed only by ∆Φ = AZ av . Theprediction Z c ( θ M ) = AZ av from the theory (Eq. (13))fits well with the simulation results. In contrast, with astronger impulse (Fig. 3(c)), the simulation results de-viate from Z c ( θ M ) = AZ av for some range of θ M . Wecalculate h ∆ φ i from the simulations and it accounts forthe deviation as predicted from Eq. (9).The deviation in Fig. 3(c) can be understood as fol-lows. Nonzero ∆ φ i can occur only when the order pa-rameter transiently decreases. For θ M ∈ ( π/ , π ), theimpulse disperses the locked group, because the trailingoscillators receive more negative impact than the lead-ing ones. Thus, the order parameter decreases from R to R ( t +) ( < R ). R can decrease more depending onthe behavior of the oscillators, and then returns to R .The collective phase Θ also decreases ( ˜Θ( t +) <
0 400 800 1200 1600 2000 (b)(c) θ π M =0.825-0.5-0.4-0.3-0.2-0.1 0 0.1 0.2 0.3 0.4 0.5 ZAZAZAZ − ∆ φ < > c1av av sim
R(t +) - R -0.1 0 0.1 0.2 0.3 0.4 0.5 -0.5-0.4-0.3-0.2-0.1 0 0.1 0.2 0.3 0.4 0.5 -1.5-1-0.5 0 0.5 1 1.5 θ M θ M θ M θ M ∆ φ i θ i (a)(c) (b)(d)(e) (f ) PR C PR C Osc Index i Osc Index i PR C PR C FIG. 3: (Color online) Fully locked cases. (a) snapshot ofphases of oscillators: a fully locked state (b) A = 0 .
5. (c) A = 1 .
0. (d) ∆ φ i for (b) and (c) with θ M = 0 . π . For(a)-(d), H ( θ ) = sin θ and Z ( θ ) = − sin θ . (e) H ( θ ) =sin θ , Z ( θ ) = − [sin( θ + 0 . π ) − sin(0 . π )] / [1 + sin(0 . π )],and A = 0 .
5. (f) H ( θ ) = sin( θ + 0 . π ), Z ( θ ) = − sin θ ,and A = 0 .
5. A uniform distribution is used for { ω i } : (i) ω i = ¯ ω − ∆ ω + ω ( i − N − or (ii) randomly selected ω i from[¯ ω − ∆ ω, ¯ ω + ∆ ω ]. (i) is used for the symbols of (a)-(e) exceptthe gray circles in (b), (c), (e), and (f). ¯ ω = π/
2. ∆ ω = 0 . K = 0 . K = 1 for (f). to the impulse. The sudden changes in R and ˜Θ af-fect the dynamics of oscillators. The behaviors of oscil-lators right after the impulse can be described by theequation ˙ ψ i = ω i − Ω + KR ( t +) sin( ˜Θ( t +) − ψ i ) with ψ i ( t +) = ψ i ∗ + Z ( θ M ( t ) + ψ i ∗ − ψ M ∗ ), where ψ i ∗ is thelocking phase for oscillator i . The trajectory of oscillator i can escape completely from the basin of attraction of φ i during the transient behavior of R and settle to theequivalent phases φ i = φ i + 2 n i π . Since the curves for( ψ i , ˙ ψ i ) are shifted to the left due to the negative ˜Θ( t +)and upwards(downwards) for the oscillators with ω i > Ω( ω i < Ω), the oscillators with frequencies far from theaverage one can escape and those with higher frequen-cies escape first. Because of this, h ∆ φ i > AZ av (Eq. (13)). The oscilla-tors with frequencies far from the average one have morechance to have higher n (Fig. 3(d)), because they candrift faster and stay unlocked longer. Other ranges of θ M can be understood similarly.When H ( θ ) is not odd, it is difficult to find the gen-eral results. For H ( θ ) = sin( θ + β ), we can see that thesecond term of Eq. (10) is equal to KR sin β . Thus,∆Φ = AZ av + K sin β R ∞ t ( R − R ) dt . With weak stimulus, Z c ( θ M ) = ∆Φ and the PRC deviates posi-tively (negatively) from AZ av for values of R ( t +) > R =0 M =0.4 M =0.55 M =0.85 M −0.6−0.4−0.2 0 0.2 0.4 0.6 2π3π/2ππ/20 ZAZAZAZ − ∆ φ < c1av av sim R(t +) - R > θ M θ i ∆ φ i (a)(b) (c) PR C Osc Index i Osc Index i
FIG. 4: (Color online) Partially locked cases with H ( θ ) =sin θ . (a) snapshot of phases of oscillators: a partially lockedstate. (b) A = 0 .
5. (c) ∆ φ i . Z ( θ ) = − sin θ . { ω i } obeys a Gaussian distribution g ( ω, ¯ ω ) = σ √ π exp( − ( ω − ¯ ω ) σ ):(i) ω ( N/ k = ¯ ω + y k , ω ( N/ − k +1 = ¯ ω − y k with y k =( x k − + x k ) / x k +1 = x k + N − /g ( x k , x = 0 for k = 1 , ..., N/ ω i = ¯ ω − y i and ω i + N/ = ¯ ω + y i with i ≤ N/ y i randomly selected according to g ( y, y >
0. (i) is used for the simulations(symbols) of (a)-(c)except the gray circles in (b). ¯ ω = π/ σ = 0 .
3, and K = 0 . ( R ( t +) < R ) (Fig. 3(f)). The values of R ( t +) areeasily calculable using Z and the distribution for ω .Finally, let us briefly consider partially locked cases inthe limit of N → ∞ . When the system exhibits a par-tially locked state, the drifting oscillators form a station-ary distribution in the frame rotating with the synchro-nization Ω. In the original frame, we can say that the dis-tribution rotates with Ω. We can define PRCs for lockedoscillators. Let us consider cases with H ( θ ) = sin θ . Fig-ure 4(b) shows the PRC Z c ( θ M ) of oscillator M for thepartially locked state of (a). We can understand Z c ( θ M )through Eq. (9). Since H is an odd function, ∆Φ = AZ av (Eq. (13)). While for the locked oscillators ∆ φ i = 2 n i π and is nonzero in some ranges as in the fully locked cases,for the drifting oscillators ∆ φ i is usually not an integermultiple of 2 π and can be nonzero in any ranges (Fig.4(c)). Simulations show that ∆ φ i for the drifting oscil-lators contribute significantly to the PRC and the PRC(symbols, Fig. 4(b)) differ from AZ av (the dashed curve)for almost the entire range of θ M .In summary, we have investigated the PRCs of coupledoscillators in terms of the PRCs of individuals, the natureof the coupling, and the relative phases of the oscillators.Our approach of obtaining PRCs using the average andrelative phases can be applicable to oscillators on differ-ent type of networks.This work was supported by National Science Founda-tion grant DMS05135. [1] A. Pikovsky, M. Rosenblum, and J. Kurths, Synchroniza-tion: a universal concept in nonlinear sciences (Cam-bridge University Press, Cambridge, 2001)[2] Y. Kuramoto,
Chemical Oscillations, Waves, and Turbu-lence (Springer, Berlin, 1984); S. H. Strogatz, Physica D , 1 (2000); J. A. Acebr´on et al. , Rev. Mod. Phys. ,137 (2005).[3] A. T. Winfree, The Geometry of Biological Time , 2nd ed.(Springer, New York, 2001).[4] J. Rinzel and G. B. Ermentrout, in
Methods in NeuronalModeling , 2nd ed. (MIT Press, Cambridge, MA, 1998).[5] G. B. Ermentrout and D. Kleinfeld, Neuron , 33(2001).[6] P. Tass, Phase Resetting in Medicine and Biology (Springer, Berlin, 1999).[7] G. B. Ermentrout, Neural Comp. , 979 (1996); L. Glass,Y. Nagai, K. Hall, M. Talajic, and S. Nattel, Phys. Rev. E , 021908 (2002); E. Brown, J. Moehlis, and P. Holmes,Neural Comp. , 673 (2004); R. F. Gal´an, G. B. Er-mentrout, and N. N. Urban, Phys. Rev. Lett. , 158101(2005); R. Gunawan and F. J. Doyle III, Biophys. J. , 2131 (2006); E. M. Izhikevich, Dynamical Systemsin Neuroscience (MIT Press, Cambridge, MA, 2007).[8] Y. Kawamura et al. , Phys. Rev. Lett. , 024101 (2008).[9] H. Sakaguchi and Y. Kuramoto, Prog. Theor. Phys. ,576 (1986).[10] H. Daido, Phys. Rev. E61