Spatial Effects of Phase Dynamics on Oscillators Close to Bifurcation

: The phase reduction approach has manifested its efficacy in investigating synchronization behaviors in limit-cycle oscillators. However, spatial distributions of the phase value on the limit cycle may lead to illusions of synchronizations for oscillators close to bifurcations. In this paper, we compared the phase sensitivity function in the spatial domain and time domain for oscillators close to saddle-node homoclinic (SNH) bifurcation, also known as saddle-node bifurcation on an invariant circle. It was found that the phase sensitivity function in the spatial domain can show the phase accumulation feature on the limit cycle, which can be ignored in the phase sensitivity function in the time domain. As an example, the synchronization distributions of uncoupled SNH oscillators driven by common and independent noises were investigated, where the space-dependent coupling function was considered on common noise. These results shed some light on the phase dynamics of oscillators close to bifurcations.


Introduction
The phase reduction approach [1] is an efficient method to investigate the dynamical behaviors of periodic oscillators subjected to weak perturbations by utilizing the phase variable to reduce the dimensionality of the original system.As the key part, the phase sensitivity function (PSF) [2] or the infinitesimal phase response curve (PRC) [3] directly reflects the impact of the perturbation on the change of the phase value (advance or delay).The PSF has a deep connection with the conception of isochrons [4,5] and Floquet eigenvectors [2].Due to its simplicity and efficiency, this approach has been widely employed in coupled and uncoupled oscillators, delayed systems, quantum synchronization, network dynamics, reaction-diffusion systems, hybrid systems, relaxation oscillators, noisy limit-cycling systems, and even noise-induced coherent excitable systems [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21], etc.In stochastic dynamics, a similar idea termed the stochastic sensitivity function is proposed [22], which can be applied to both stable equilibrium points and limit cycles [23,24].The difference is that the phase sensitivity function considers the phase response against any perturbation, whether random or deterministic, and is irrelevant to the perturbation strength, while the stochastic sensitivity function focuses on random perturbations and depends on the noise intensity and a predefined confidence probability.
In neuroscience, Hansel et al. [25] first found two types of PRC.Type I always advances the phase via excitatory postsynaptic potentials, while Type II can either delay or advance the phase.Later, Ermentrout [26] related the types of PRC with the classification of the excitable membranes by Hodgkin and further found that Type II is more favorable for stochastic synchronization [27].Brown et al. [28] made a probabilistic analysis of phase response for four kinds of neuron models encompassing four generic bifurcations and obtained the scaling behaviors near the bifurcations, which contained both types of PRC.
As a paradigmatic model of Type I PRCs, oscillators close to saddle-node homoclinic (SNH) bifurcation have been analyzed using the phase reduction approach [3,26,28].Different from the Type II PRC, the phase of the SNH oscillator can accumulate around the bifurcation point.However, this cannot be unveiled by the traditional phase response curve or the phase sensitivity function, as the PSF is usually measured on the phase (time domain), not space.Therefore, to exhibit the spatial accumulation characteristics, we will compare the PSF both in the time domain and in the spatial domain.
This paper is organized as follows.In Section 2, the mathematical model and its bifurcation behavior will be discussed.The PSF in the time domain and spatial domain will be compared in detail in Section 3. In Section 4, the noise-induced synchronization of uncoupled SNH oscillators will be investigated for constant and space-dependent coupling functions on common noise, where the spatial effects of the couplings will be considered on the synchronization behaviors.Finally, the discussions and conclusions are given in Section 5.

Mathematical Model and Setup
We consider the canonical oscillator close to the saddle-node homoclinic bifurcation given by [29] Under the polar transformation  =  cos() ,  =  sin(), system (1) can be transformed as follows: ̇= (1 −  2 ), where r represents the radius and  is the angle variable (do not confuse this with the phase value defined in the reduced phase equations). is the bifurcation parameter.The system will always have an unstable equilibrium (0, 0) regardless of the value of .For −2 <  < 0 , the system will have two other equilibrium points: r = 1 and  = ± cos −1 (−(1 + )) (black and white dots in Figure 1).For  = 0, system (1) will have a saddle-node point ( = 1,  = ) and undergo saddle-node homoclinic bifurcation.For  > 0, a stable limit cycle appears.It should be noted that the circle with radius r = 1 is an invariant cycle for all values of .The bifurcation scenario is illustrated in Figure 1 for three typical values of .

Figure 1.
Bifurcation scenario for the saddle-node homoclinic bifurcation of system (1).The black dot represents the stable equilibrium and the white dots denote the unstable ones.The gray dot for  = 0 is the saddle-node point.The blue and red curves are the x-and y-nullclines, respectively.The gray arrows denote the directions of the vector field, and black curves display trajectories of the system with the black arrows showing the directions of the motions.
To have a stable limit cycle, we set  > 0 in the following and investigate cases when  is close to the bifurcation value (i.e.,  = 0).For simplicity, we named the oscillator close to the saddle-node homoclinic bifurcation as the SNH oscillator.It can be checked that the x-nullcline is always pinned at (±1, 0) and the y-nullcline is always pinned at (0, ±1).This property will be useful later for determining the positive or negative phase response at these points.The angle frequency of the limit cycle is  = √( + 1), and the period is 2/.It can be seen that as  → 0, the period will approach infinity.The timeangle curve can be also obtained from system (2) as The time-angle curves for system (2) are illustrated in Figure 2 for varying bifurcation parameters  .The period of the oscillator increases as  → 0 and the major time course of the trajectory remains close to angle  = , which is consistent with the saddlenode point when  = 0.

Phase Sensitivity Function in the Time Domain
Considering the weakly perturbed situation, to simplify the notation, we rewrite system (1) as where  = [, ]  , () = [  (),   ()]  is the original vector field, and (, ) is the weak perturbation.Using the standard phase reduction approach [1,30], it can be reduced to the system with a single phase variable as ⋅ ( 0 (()), ) =  + () ⋅ (, ), where  0 () is the limit cycle of the unperturbed system and () = [  (),   ()] is the socalled PSF or the infinitesimal PRC, which is the gradient of the phase variable on the limit cycle [1][2][3].The phase dynamics of system (1) can be captured through the reduced phase equation.
Although PSF cannot be analytically calculated in most cases, there are a handful of numerical methods.Here, we apply two methods, namely, the direct method and the adjoint method [1,3,26], to calculate the PSF of the SNH oscillator.Figure 3 illustrates the results for  = 0.2 using these two methods, which are consistent with each other.Phase sensitivity function (time domain) for the SNH oscillator, calculated using the direct method (red dots) and the adjoint method (black line), respectively.The parameter  = 0.2.In the direct method, the perturbation intensity was fixed at 0.005.The initial phase  = 0 was chosen at (0, 1).(a) X-component; (b) Y-component.
For  → 0, the results for PSF have been investigated widely [3,26,28].However, for comparison, we still calculated them using the adjoint method as in Figure 4.By decreasing , the PSF gradually approaches the curve  sin 2 () or (1 − cos(2)) [3,26], where A and K are proper constants.The limit situation shows that the PSF along the x-direction, i.e.,   (), will be nonnegative (similarly, nonpositive for   ()).This is called the Type I PRC by Ermentrout [26].However, it should be noted that for a smooth two-dimensional dynamical system, it is impossible for the PSF or PRC to be strictly nonnegative.A nonnegative PSF is only possible for the limit situation, e.g., by making  = 0 in the SNH oscillator (but this will cause the period of the limit cycle to be infinity) or for a non-smooth or discontinuous system (e.g., the quadratic integrate-and-fire neuron [3]).This property can be verified using the phase sensitivity function in the spatial domain discussed in the next subsection.

Phase Sensitivity Function in the Spatial Domain
The phase sensitivity function in the time domain (PSFt) for the SNH oscillator showed an overall increase in amplitude when the parameter was close to the bifurcation (see Figure 4).Additionally, the nonnegative interval for   () extended gradually to the whole period.However, at the same time, as far as the space was concerned, the phase was compressed around the angle  = .Figure 5 shows this feature by depicting the phase values on the limit cycle.Although the phase by definition uniformly increases with time, it is not uniform in space.As  → 0, the phase accumulates around the phase  =  or the space position (−1, 0).To illustrate the phase sensitivity function in the spatial domain (PSFs), we reparametrized the phase sensitivity along the limit cycle as a function of the space position.Similar to the PSFt, we may define a space variable  which varies from 0 to 2.By assumption,  will be uniform along the limit cycle within which  will increase linearly as the arc length (i.e., space).Because the limit cycle in our case is just a unit circle centered at the origin, by Equation (3) (note that the starting angle of Equation ( 3) is 0, but the starting angle in Figure 5 is /2), the phase sensitivity in the spatial domain ( in our case) can be obtained.The PSFs is plotted in Figure 6, where the parameter  is chosen the same as that in Figure 4.By comparing the PSFs with the PSFt, it is noted that the phase sensitivity is broad in time but narrow in space.The PSFs is peaked mainly around the angle  = , which is more significant as  decreases.This is consistent with the results in Figure 5, where the phase value accumulates at  = , so the phase change will be more sensitive to the perturbation there.Next, we want to obtain some quantitative and qualitative characteristics of the PSFs.Considering the definition of the phase function and the phase sensitivity function, the following condition always satisfies: which is also the normalization condition for the adjoint equation [1,26].It is the reason for the origin (0, 0) to be the phaseless equilibrium [3,5], as the vector field vanishes there.
For the SNH oscillator, as is mentioned previously, there are four intersection points for the limit cycle and the nullclines: (±1, 0) for the x-nullcline and (0, ±1) for the y-nullcline.
In general, for any smooth two-dimensional dynamical system, and for a perturbation in any direction, say  ⃗ , the limit cycle will have at least two positions.For one position, the vector field on it will be along  ⃗ and the other will be against  ⃗ (or along the direction − ⃗ ).Thus, using Equation ( 5), the PSFs along  ⃗ will be positive at the former position while negative at the latter.Therefore, for any smooth two-dimensional dynamical system, strictly nonnegative PSF (or Type I PRC) is impossible; there are at least two positions having PSF with opposite signs (and their neighboring area by continuity), although the amplitude may be small.In other words, every smooth two-dimensional limitcycling system possesses Type II PRC. Figure 7 illustrates a schematic diagram for the explanation.For systems with higher dimensions, it remains to be validated whether there is a similar result as in the two-dimension case.In the following section, we will show that PSFs can reveal the importance of spacedependent perturbation on the system, while PSFt is important to visually show the degree of synchronization.

Spatial Effects of Synchronization Using Space-Dependent Coupling
In this section, we consider uncoupled SNH oscillators driven by common and independent noises.The governing equation is as follows: where () is the common noise shared by all the SNH oscillators and   () is the independent noise.They are assumed to be independent, identically distributed zero-mean Gaussian white noise and satisfy the correlation functions ⟨  ()  ()⟩ =  , ( − ) , ⟨   ()   ()⟩ =  ,  , ( − ), and ⟨  ()   ()⟩ = 0. D1 and D2 denote their strengths.According to [7], derived by Nakao et al., the stationary probability density distribution (PDF) of the phase difference  =   −   is given as , where u0 is the normalization constant.The functions () and () are given as () = ()  ( + )  ( + ).Additionally, (0) and (0) are corresponding values at  = 0.In the independent noise case, i.e.,  1 = 0,  2 ≠ 0 , the PDF of the phase difference will be a uniform e e − distribution: () = 1/2, whereas under the common noise condition, i.e.,  1 ≠ 0,  2 = 0, the PDF will be a delta function with zero phase difference.This is the phenomenon of common noise-induced complete synchronization.In this paper, we consider the existence of both kinds of noises.
We set noise strengths as  1 = 0.005,  2 = 0.0001.Figure 8 illustrates the PDF of the phase difference and by Monte Carlo simulations for different coupling functions ( () ) on the common noise.For Figure 8a, ( () ) =  1 = diag(1, 1), and the common noise is uniform in space.The theoretical and numerical results show that the PDF of the phase difference is peaked at  = 0, which is a typical one-cluster synchronization.The inset shown in the top left corner of Figure 8a illustrates the spatial distribution of the oscillators on the limit cycle.For space-dependent coupling functions, we analyze three kinds of exponential kernels, which are: (1) ( () ) =  2 =  −((+1) 2 + 2 ) ⋅ diag(1, 1) ; ( 2) ( () ) =  3 =  −( 2 + 2 ) ⋅ diag(1, 1) ; and (3) ( () ) =  4 =  −((−1) 2 + 2 ) ⋅ diag(1, 1) .Figure 8b-d show the PDF for three different kernels on the common noise, peaking at (−1, 0), (0, 0), and (1, 0), respectively.The PDFs for Figure 8a,b are almost the same.This shows that the constant coupling function ( () ) and the one with the exponential kernel centered at (−1, 0) give rise to similar synchronization behaviors for the SNH oscillators.However, for exponential kernels centered far from (−1, 0), the PDF can be much flatter, which shows a rather weaker synchronization (Figure 8c), and can even induce asynchrony with a nearly uniform distribution (Figure 8d).The above observation can be simply inferred from the PSFs.As we can see from the PSFs in Figure 6, the PSFs is peaked around  =  or (−1, 0).Therefore, the perturbation around this position can have a more significant influence on the phase response and the synchronization behaviors.
For other values of the bifurcation parameters , we calculate the PDFs in Figure 9 for the same four kinds of coupling functions on the common noise.The results show a similar tendency (only a slight difference is observed, such as the insets in Figure 9).It should be noted that although the instantaneous distributions of the SNH oscillators on the limit cycle can display, to some extent, the synchronization behaviors of the ensembles, they may lead to some illusions only based on them.Figure 10 compares the snapshots of the instantaneous distribution of oscillators on the limit cycle (spatial domain) and on the phase circle (time domain) for different parameters and coupling functions.From the snapshots for the oscillators on the limit cycle, at each value of , the degree of synchronization gradually attenuates from left to right as expected.However, for the same coupling function, the instantaneous distribution of the SNH oscillators exhibits a more "clustered" behavior for a smaller value of  (the third column in Figure 10), which does not imply a more synchronized state.On the contrary, the PDF demonstrates that as  decreases, the peaks reduce slightly (see Figure 9), which leads to an opposite conclusion.It can be explained by the spatial nonuniform distribution of the phase value, as is previously revealed in Figure 5.This is quite different from the uniform situation, where snapshots of the instantaneous distribution of the oscillators can clearly perform their clustering features, as in Ref. [7] for the Stuart-Landau oscillators.On the other hand, as  → 0, the snapshots of the oscillator distribution on the limit cycle will give less information compared with the stationary PDF (see, e.g.,  = 0.01, where the oscillators are "clustered" even when the PDF is almost uniform).The above results show that for oscillators such as the SNH system close to bifurcation, the accumulated phase value on the limit cycle may mislead us on the systems' synchronization behaviors.In that case, the phase circle where the phase value is uniformly distributed can display the true degrees of synchronization (see Figure 10 for the oscillators on the phase circle).

Conclusions
In summary, the phase sensitivity functions in the time domain and spatial domain of oscillators close to saddle-node homoclinic bifurcation are investigated in this paper.It is found that for the SNH oscillator close to bifurcation, the phase sensitivity function in the time domain (PSFt) is close to the quadratic sine curve, which is a well-known result as the Type I PRC.However, the phase sensitivity function in the spatial domain (PSFs) indicates that there will always be positions on the limit cycle for PSF to be positive and negative.A nonnegative PRC can only be realized in the limit situation or for non-smooth or discontinuous dynamical systems.That means that for any smooth two-dimensional limit-cycling system, Type II PRC is the only possible type.Additionally, the PSFs can display the phase accumulation effect on the limit cycle, which can explain the different synchronization behaviors by space-dependent perturbations.It should be noted that, due to the phase accumulation effect, the population distribution on the limit cycle may give illusions on synchronizations.In that case, it is suggested that the phase circle should be used to display the real synchronization characteristics.This is consistent with a recent work by Freitas et al. [31], wherein they found that topological properties of an attractor have no bearing on phase coherence.
In future works, spatial effects on the phase dynamics can be investigated in higher dimensional systems with more complex bifurcation structures, wherein multiple timescales may play a vital role.Additionally, large-scale networks with complex connection topologies can imply more possibilities between spatial and temporal effects.

Figure 2 .
Figure 2. Time-angle curves for one period of system (2) with different bifurcation parameters .
, ), where () = Θ(()) is the phase variable of the system, which by definition linearly grows with time t as () = .Note here that the phase variable () is different by definition from the angle variable () defined previously.They are linear in time and in space, respectively.Due to the weak perturbations, it can be further reduced to the linear order of the perturbation as

Figure 3 .
Figure 3.Phase sensitivity function (time domain) for the SNH oscillator, calculated using the direct method (red dots) and the adjoint method (black line), respectively.The parameter  = 0.2.In the direct method, the perturbation intensity was fixed at 0.005.The initial phase  = 0 was chosen at (0, 1).(a) X-component; (b) Y-component.

Figure 4 .
Figure 4. Phase sensitivity function (time domain) for the SNH oscillator for varying bifurcation parameters .The solid curves and the dashed curves represent results for   () and   (), respectively.The red and green arrows show the change in the PSF by decreasing .The blue dashdotted line is the zero-value line.Parameters are:  = 0.2, 0.1, 0.07, 0.05, 0.03, 0.01.

Figure 5 .
Figure 5. Phase value  on the limit cycle for different bifurcation parameters.The line segments are phases equally distributed along the limit cycle with the interval /20.

Figure 6 .
Figure 6.Phase sensitivity function (spatial domain) for the SNH oscillator for varying bifurcation parameters .(a) The solid curves and the dashed curves represent results for   () and   (), respectively.(b) The three-dimensional view for PSFs.Parameters are the same as those in Figure 4.

Figure 7 .
Figure 7. Schematic diagram of positions for phase advance and delay with perturbation direction  ⃗ .

Figure 9 .
Figure 9. Stationary probability density function of the phase difference for different bifurcation parameters.The four kinds of coupling functions are the same as those in Figure 8.The insets display the local enlargements.

Figure 10 .
Figure 10.Snapshots of instantaneous distribution of oscillators on the limit cycle and on the phase circle for different parameters and coupling functions.