Evaluating the phase dynamics of coupled oscillators via time-variant topological features
EEvaluating the phase dynamics of coupled oscillatorsvia time-variant topological features
Kazuha Itabashi, ∗ Quoc Hoan Tran, † and Yoshihiko Hasegawa ‡ Graduate School of Information Science and Technology,The University of Tokyo, Tokyo 113-8656, Japan (Dated: May 11, 2020)The characterization of phase dynamics in coupled oscillators offers insights into fundamentalphenomena in complex systems. To describe the collective dynamics in the oscillatory system, orderparameters are often used but are insufficient for identifying more specific behaviors. We thereforepropose a topological approach that constructs quantitative features describing the phase evolutionof oscillators. Here, the phase data are mapped into a high-dimensional space at each time point, andtopological features describing the shape of the data are subsequently extracted from the mappedpoints. We extend these features to time-variant topological features by considering the evolutiontime, which serves as an additional dimension in the topological-feature space. The resulting time-variant features provide crucial insights into the time evolution of phase dynamics. We combinethese features with the machine learning kernel method to characterize the multicluster synchronizeddynamics at a very early stage of the evolution. Furthermore, we demonstrate the usefulness of ourmethod for qualitatively explaining chimera states, which are states of stably coexisting coherentand incoherent groups in systems of identical phase oscillators. The experimental results show thatour method is generally better than those using order parameters, especially if only data on theearly-stage dynamics are available.
I. INTRODUCTION
Coupled phase oscillators have been widely used to in-vestigate cooperative behaviors in complex systems. Thecoupling scheme of the oscillators in nature are reflectedby the dynamic behaviors of complex systems, whichinclude chimera state emergence, chaos, multistability,and synchronization. For example, a brain network canbe analyzed as a coupling model of ten billion neurons,which can exhibit chimera-like states with responses re-lating to certain brain disorders [1–3]. Meanwhile, circa-dian clock systems, in which a synchronized state and anasynchronized state interchange within a one-day cycle,are related to other diseases such as metabolic disorders,cataplexy, and narcolepsy [4, 5]. Electromechanical sys-tems can also behave as chaotic oscillators, and gainingan appropriate understanding of oscillator dynamics thuspromotes the realization of these systems, which can op-erate in a noisy environment [6, 7]. One of the primarypoints of interest related to coupled oscillators is the pre-diction of the attendant dynamics from data obtainedat the early stage of the system. Because many collec-tive dynamics can be modeled by coupled oscillators, theearly prediction of coupled oscillator dynamics has goodpotential for applications such as the diagnosis of humandiseases and the detection of specific malfunctions. How-ever, the realization of these applications requires theadoption of theoretical and computational methods toadequately represent the time profile of coupled oscilla-tors. ∗ [email protected] † [email protected] ‡ [email protected] The intrinsic dynamic properties of coupled oscillatorsare often described in terms of their phase variables.Meanwhile, the degree of the global synchronization ofthe system is commonly expressed by the Kuramoto or-der parameter, which represents the phase coherence ofthe oscillators. The Kuramoto order parameter takes thevalue of 0 for complete asynchronization and up to 1 forfull synchronization [8]. However, the order parametermethod is insufficient for effectively analyzing the syn-chronization situation under certain conditions. For ex-ample, when the value of the order parameter is 0, theoscillators can be synchronized in terms of a symmetricalphase distribution even though they should, by defini-tion, be completely asynchronous. In view of this, ratherthan use specific global parameters to represent the col-lective dynamics, our approach focuses on the topologi-cal aspects of specific phase variables to reveal the phasedynamics along the evolution timeline in terms of moreefficient quantities. First, at each oscillator, we define aphase-dependent point in a high-dimensional space. Wethen focus on the time-variant evolutionary change inthe shape of these points known as a point cloud P ( t )within the space. We hypothesized that this evolutionarychange has a close relationship with pattern formation,signal propagation, and the stochastic phenomena andextensive chaos in oscillatory systems. Thus, we demon-strate that tracking these changes can provide crucial in-sights into the dynamics at the early stage of the system.The basic idea of our approach is that topological as-pects can help to reveal the underlying structure of thecollective patterns of oscillators as the system evolves.Here, we apply persistent homology analysis [9] to eval-uate the configuration of P ( t ) in terms of quantitativetopological features and to monitor the variation of thesefeatures throughout time t for application in the task of a r X i v : . [ phy s i c s . d a t a - a n ] M a y predicting the dynamics. Persistent homology is a tech-nique from algebraic topology in approaches aimed atrepresenting the shape of the data related to the infor-mation on the topological structures, such as the con-nected components, loops, or holes in the data. Persis-tent homology has been effectively applied to studyingthe qualitative changes in data from various dynamicalsystems, including time series [10, 11], time-varying net-works [12–14], and quantum data [15]. Given a non-negative threshold ε , we place ε -radius balls centered atpoints in P ( t ) and observe the shape of the space over-lapped by these balls [Fig. 1(a)]. At a sufficiently small ε , we can obtain the shape without any difference be-tween the original points in the point cloud. When weincrease ε , the balls intersect, which subsequently leadsto the change in the topological structures in terms of,for example, the connected components or loops in thespace. Here, the connected components tend to merge,while loops emerge and then vanish in accordance withthe gradual change in the threshold. At each time point,we define the topological features as the values of ε torepresent the emergence and disappearance of the topo-logical structures [Fig. 1(b)]. Furthermore, we extendthese features by adding the time axis [Fig. 1(c)].These extended features, referred to as time-varianttopological features, can reflect the temporal behavior ofthe oscillators and can therefore provide useful knowledgefor predicting the dynamics at the early stage of the sys-tem. In fact, these features can serve as discriminate fea-tures in qualitative evaluations of the phase dynamics ofoscillators. We can input these features into the machinelearning kernel algorithms to allow for the applicationof statistical-learning tasks, which include classifying thebehaviors of multicluster synchronization or predictingchimera states in which synchronization and asynchro-nization co-exist. Interestingly, our approach character-izes the phase dynamics of oscillators at a very early stageof the time evolution, a point where conventional orderparameters do not effectively operate. II. METHODSA. Time-variant topological features
To explain the proposed time-variant topological fea-tures, we considered the Kuramoto model, which is themost common and best-suited model for understandingsynchronization phenomena in physical, chemical, andbiological systems. Given a set of N coupled heteroge-neous oscillators, the Kuramoto model was formalized asa set of first-order differential equations: dθ i dt = ω i + 1 N N (cid:88) j =1 g ij sin( θ j − θ i − α ) . (1)Here, θ i and ω i denote the phase and the natural fre-quency of the oscillator i , respectively, and g ij ≥ i and j .The angle α is a tunable parameter describing the phaselag between oscillators i and j .To study the dynamics of the coupled oscilla-tors, we considered a mapping ϕ from the set { θ ( t ) , θ ( t ) , . . . , θ N ( t ) } of phases to the set P ϕ ( t ) = { s ( t ) , s ( t ) , . . . , s N ( t ) } of points in a specific L -dimensional space R L , where s i ( t ) = ϕ ( θ i ( t )) ∈ R L . Byconsidering an appropriate mapping ϕ , the informationon the evolutionary change in the configuration of P ϕ ( t )along t can provide important insights into the dynamicsof coupled oscillators. Our approach to quantifying theconfiguration of P ϕ ( t ) into quantitative features relies onthe persistent homology theory. First, we defined a dis-tance function d ϕ : P ϕ ( t ) × P ϕ ( t ) → R to evaluate thedissimilarity between oscillators i and j at time t . Wecentered a ε -radius ball at each point s i ( t ) in P ϕ ( t ), i.e.,to form the set B ( ε, s i ( t )) = { v ∈ R L | d ϕ ( s i ( t ) , v ) ≤ ε } .Then, by taking the union of these balls over all i =1 , , . . . , N , we could obtain an overlapped space, U ( ε, P ( t )) = N (cid:91) i =1 B ( ε, s i ( t )) , (2)the shape of which represents the configuration of P ϕ ( t )at the radius ε . The idea of persistent homology is totrack how this shape changes as the radius ε increases.In fact, the method allows for modeling the shape of U ( ε, P ϕ ( t )) in a far more mathematically and computa-tionally tractable representation, i.e., a simplicial com-plex, which is a complex of geometric structures knownas simplices. Here, an n -simplex presents a generaliza-tion of the notion of a triangle or tetrahedron to arbi-trary dimensions. For example, a 0-simplex is a point, a1-simplex is a line segment (with two end points as itsfaces), and a 2-simplex is a triangle and its enclosed area(with three edges and three vertices as its faces). Simi-larly, a 3-simplex is a filled tetrahedron (with triangles,edges, and vertices as its faces), and while a 4-simplexis beyond visualization, it is a filled shape with tetrahe-drons, triangles, edges, and vertices as its faces. One ofthe main types of simplicial complexes is known as theVietoris–Rips complex, which we will now briefly review.A Vietoris–Rips complex V R ( ε, P ϕ ( t )) is a collection ofsimplices, where each simplex is built over a subset ofpoints in P ϕ ( t ) if B ( ε, s i ( t )) ∩ B ( ε, s j ( t )) (cid:54) = ∅ for everypair of points s i ( t ) , s j ( t ) in the subset. Now, startingwith ε = 0, the complex contains only the 0-simplices,i.e., the discrete points. As ε increases, connectionsemerge between the points and the edges (1-simplices),and filled triangles (2-simplices) are introduced into thecomplexes [Fig. 1(a)]. This process enables us to obtaina sequence of embedded complexes, which is defined as afiltration. Moreover, if ε becomes considerably large, allthe points become interconnected, which means no usefulinformation can be conveyed.As noted above, persistent homology is focused on theemergence and disappearance of topological structures (a) (b) D ea t h Birth ε = 0 ε = 1 ε = 2 ε = 3 ε = 4(c) π-π0 ………… τ k τ i τ j τ k τ i τ j T i m e B i r t h D ea t h t t Connected componentsLoops
Persistence diagram D (3) Mapping φθ i ( t )→ s i ( t )Time series data φ ( τ k ) φ ( τ j ) φ ( τ i ) ε Time θ i ( t )- θ ( t ) FIG. 1. Illustration of time-variant topological features for coupled oscillator systems. (a) Example of a sequence of theVietoris–Rips complex, i.e., a filtration constructed from a set of discrete points. Each ball of radius ε was placed at each point,then, the Vietoris–Rips complex was used to model the shape of the union of these balls. In the complex, each simplex is formedfrom a subset of points if every pair of corresponding balls intersect. The evolution changes in the topological structures, e.g.,the merging of connected components and the emergence and disappearance of loops, were tracked through increasing the ε until no change was observed. At ε = 0, there were nine points corresponding with nine connected components in the space.At ε = 2, several edges were added into the complex, at which point a number of connected components perished and merged.Similarly, two loops emerged in the space at birth radii ε = 2 and ε = 3, respectively, and then perished at death radii ε = 3 and ε = 4, respectively. (b) The collection of birth and death radii is represented in a two-dimensional persistence diagram . (c) Ateach time step, the coupled oscillator phases were mapped into a set of discrete points, i.e., a point cloud within a metric space.The two-dimensional persistence diagram was constructed from this point cloud for each time step. Time-variant topologicalfeatures were obtained by concatenating these diagrams across the time steps into a three-dimensional persistence diagram. such as any connected components and loops in the fil-tration. For example, in Fig. 1(a), while there are nineconnected components in the space at ε = 0 and ε = 1, at ε = 2, several components are merged together, meaningthree connected components remain in the space. Simi-larly, one loop exists in the space from ε = 2 to ε = 3,and one in the space from ε = 3 to ε = 4. These topo-logical structures are mathematically represented by theconcept of the dimension 0 and dimension 1 persistenthomology groups, which are vector spaces with dimen-sions corresponding to the number of connected compo-nents and the number of loops, respectively (see [16] formore details). We use the emergence and disappearanceof these topological structures in the filtration to quan-tify the evolution in the shape of U ( ε, P ϕ ( t )). Here, eachof the connected components and loops was assigned toa persistence pair of radius ( ε b , ε d ), such that they origi-nated at the birth radius ε = ε b and perished at the deathradius ε = ε d . The collection of all persistence pairsin a two-dimensional coordinate system presents a two-dimensional persistence diagram D (2) ( P ϕ ( t )), which con-tains topological features for P ϕ ( t ) [Fig. 1(b)]. To map the dynamic properties of coupled oscillators into topo-logical features, we proposed using time-variant topolog-ical features containing time-related information on thedynamics in addition to the birth radius and death ra-dius. We extended two-dimensional persistence diagramsinto three-dimensional persistence diagrams as D (3) ( ϕ ) = { ( ε b , ε d , τ ) | ( ε b , ε d ) ∈ D (2) ( P ϕ ( τ )) ,τ = τ , τ , . . . , τ T } . (3)Here, a three-dimensional persistent diagram is formedby concatenating two-dimensional diagrams along thetime-axis at time steps τ < τ < . . . < τ T [Fig. 1(c)]. B. Kernel method for topological features
In machine learning tasks using time-variant topolog-ical features, we are typically provided with a collectionof inputs D = { D , . . . , D M } from a certain set of di-agrams, from which we must quantify any patterns toevaluate previously unseen data. Since a persistence di-agram is a multiset of points of variable size, it can bedifficult to apply the algorithms when considering the in-puts as vectors or when requiring the inner product formof calculation for the input data. Here, we employedkernel methods that use a similarity measure κ ( D i , D j )between any two diagrams D i , D j ∈ D . More precisely, afunction κ : D ×D → R is called a kernel if the Gram ma-trix K with entries K ij = κ ( D i , D j ) is positive semidef-inite. To define the kernel, a feature map Φ was consid-ered in view of mapping a diagram D i ∈ D to a vectorΦ( D i ) in a Hilbert space H b wherein we could definethe inner product (cid:104)· , ·(cid:105) H b on H b . Every feature map Φdefines the kernel κ ( D i , D j ) := (cid:104) Φ( D i ) , Φ( D j ) (cid:105) H b . Sev-eral proposals have been forwarded to define a kernelfor two-dimensional diagrams, including the persistencescale space kernel [17], which is based on the heat diffu-sion kernel, the persistence weighted Gaussian kernel [18],which emerged from kernel mean embedding, the slicedWasserstein kernel under Wasserstein geometry [19], andthe persistence Fisher kernel [20], which relies on Fisherinformation geometry. Here, we extended the persistenceFisher kernel to define a kernel for three-dimensional per-sistence diagrams, which we briefly review below.The persistence Fisher kernel considers each persis-tence diagram as the sum of normal distributions andthen measures the similarity between the distributionsvia the Fisher information metric. Let D i ∆ and D j ∆ be the point sets obtained by projecting two persistencediagrams D i and D j on the diagonal line y = x of theCartesian plane. The kernel compares two extended per-sistence diagrams, D (cid:48) i = D i ∪ D j ∆ and D (cid:48) j = D j ∪ D i ∆ ,which have the same number of points. The summa-tion of normal distributions on D (cid:48) i can be defined as ρ D (cid:48) i = (cid:80) u ∈ D (cid:48) i Z N ( u , ν I ) . Here, N is a Gaussian func-tion with bandwidth ν , I is an identity matrix, and Z = (cid:82) Ω (cid:80) u ∈ D (cid:48) i N ( x ; u , ν I ) d x is the normalization constantwith the integral calculated on Ω = D i ∪ D i ∆ ∪ D j ∪ D j ∆ .Given a positive scalar α , the persistence Fisher kernel isdefined as κ F ( D i , D j ) = exp( − αd F ( D i , D j )) , (4)where d F ( D i , D j ) = arccos (cid:16)(cid:82) Ω (cid:113) ρ D (cid:48) i ( x ) ρ D (cid:48) j ( x ) d x (cid:17) isthe Fisher information metric between ρ D (cid:48) i and ρ D (cid:48) j . III. RESULTSA. Multicluster synchronization
A network system of multiple coupled oscillators candemonstrate multicluster synchronization, i.e., the net-work may split into several clusters of independent syn-chronized or organized behavior rather than form an en-tire system of synchronized behavior. Multicluster syn-chronizations are found in asymptotic states [21–23], intransient states [24], and in modular and hierarchical structures [25]. Here, we demonstrate that our time-variant topological features obtained at the early stageof the dynamics can help to predict the multicluster syn-chronized behavior of oscillators.We consider three oscillator networks with differentconfigurations of coupling strength g ij as presented inEq. (1): the globally coupled network where all couplingstrengths are equal [Fig. 2(a)], the modular coupled net-works with two [Fig. 2(b)] and four modules [Fig. 2(c)].We set the coupling strength g ij = 2 for ∀ i (cid:54) = j for theglobally coupled network. In the modular coupled net-works, the coupling strengths of the oscillators belongingto the same module were different from those belong-ing to the different modules. We set g ij = 2 for theoscillators in the same module, and g ij = 0 .
01 for thosebelonging to the different modules. Here, we set the num-ber of oscillators as N = 128, the angular frequency of i th oscillator as ω i = 1, and the tunable parameter as α = 0. Corresponding with these configurations, differ-ent behaviors of synchronization, including single-cluster[Fig. 2(d)], two-cluster [Fig. 2(e)], and four-cluster syn-chronization [Fig. 2(f)], appear as t → ∞ .We numerically solved Eq. (1) with randomly initial-ized phases θ j (0) ∈ [0 , π ) and recorded the phases θ j ( t )at each time interval ∆ τ = 0 .
8. We then used time se-quences of T T = { τ , τ , · · · , τ T } for the persistence dia-gram calculations, where τ = 0 , τ k +1 − τ k = ∆ τ ( k =1 , , . . . , T − T is the number of time steps. Tocalculate the time-variant topological features, we con-sidered the mapping ϕ : θ → (cos θ, sin θ ) to transformthe set of oscillator phases { θ ( t ) , θ ( t ) , . . . , θ N ( t ) } tothe point cloud P ϕ ( t ) = { s ( t ) , s ( t ) , . . . , s N ( t ) } , where s j ( t ) = (cos θ j ( t ) , sin θ j ( t )) lay on the unit circle in two-dimensional space. We adopted the distance functionbetween s j ( t ) and s k ( t ) as the shortest distance betweenthem measured along the unit circle.Figure 3 shows examples of the time-variant topolog-ical features obtained from the globally [Fig. 3(a)], two-module [Fig. 3(b)], and four-module coupled [Fig. 3(c)]networks. Here, the temporal transition in the dynamicscan be observed in terms of the transition in the tempo-ral patterns of the topological features. On the top rowof Fig. 3, the orange points represent the loops formedalong the phase evolution timeline. For the globally cou-pled network, the loops quickly disappeared as the os-cillators approached a synchronized state. In contrast,for the modular coupled networks, the birth radii of theloops increased as the oscillators were divided into mul-tisynchronized clusters, which were uniformly dispersed.We considered the evolution of the topological features interms of the connected components. The birth radius ofeach connected component was zero since N componentscorresponding to N oscillators appeared first. Therefore,we focused on the evolution of the death radii as illus-trated in the bottom row of Fig. 3. Here, the colored barin each plot represents the density in the distribution ofthe death radius. As represented by the right columninside each plot, at a sufficiently large time step, one t→∞t→∞t→∞ (a)(b)(c) (d)(e)(f) FIG. 2. The schematic of oscillator networks with differentcoupling configurations and the corresponding synchroniza-tion behaviors at infinite time. The vertices and edges in(a),(b), and (c) represent the oscillators and their couplingrelations. The edges with bold lines imply strong coupling( g ij = 2), while those with thin lines imply weak coupling( g ij = 0 . t → ∞ . connected component always retained the infinite valueof the death radius. To distinguish the synchronized be-havior of each coupling configuration, we examined themerging process of other connected components along theevolution timeline. In terms of the globally coupled net-work, the connected components quickly merged to formone component since the oscillators approached formingone synchronized-state cluster. At t >
8, only one com-ponent existed [Fig. 3(a)]. If the oscillators were dis-tributed in terms of multiclusters in the synchronizedstate, more connected components existed for a longerperiod of time. The number of such components corre-sponded with the number of clusters in the synchronizedstate. For example, two connected components survivedfor a period of time of t = 20 to t = 50 [Fig. 3(b)],while four connected components survived for a period of t = 30 to t = 50 [Fig. 3(c)]. These values corre-sponded with the behaviors of two-cluster [Fig. 3(b)] andfour-cluster [Fig. 3(c)] synchronization, respectively.Next, we applied the kernel method to characterizethe differences among the synchronized behaviors of os-cillators. It should be noted that our approach does notrequire prior labeling for the synchronization behaviors;rather, the focus is on characterizing the differences inthe kernel space. This approach aligns with the un-supervised learning schemes, which are fundamentallydifferent to the supervised learning schemes often usedin machine learning methods. In the supervised learn-ing schemes, the learning machine is trained on sampleswith predefined labels before the machine then attemptsto predict an unknown label of a given sample, demon-strating that it has learned by generalizing to samplesit has not encountered before. In contrast, unsupervisedschemes do not require prior labelling but characterizethe unknown dynamics via dimensional reduction meth-ods. Here, we employed the persistence Fisher kernel de-scribed in Eq. (4) for the three-dimensional persistencediagrams. The synchronized behaviors were identifiedvia the dimensional reduction method of kernel principalcomponent analysis [26].In our experiments, for each of the coupling con-figurations, we prepared 100 random initializations forthe oscillator phases. Here, different types of synchro-nization, including single-cluster, two-cluster, and four-cluster synchronization, appear as t → ∞ . We used atime sequence of T T = { τ , τ , · · · , τ T } for the persis-tence diagram calculations, where τ = 0 , τ k +1 − τ k = 0 . k = 1 , , . . . , T − T controls the periodof time to be used to detect the synchronized behaviors.We computed the Gram matrix for a total of 300 three-dimensional persistence diagrams. Figure 4(a) highlightsthe projection up to the third component of the kernelprincipal component analysis for the Gram matrix. Here,the purple, blue, and orange points represent the config-urations with single-cluster, two-cluster, and four-clustersynchronization, respectively. At the very early stage,where, for example, T = 4, it was difficult to observethe clear difference between the points belonging to dif-ferent groups of synchronized behaviors. However, theseparability increased as we increased T to T = 7 , , { r ( τ ) , r ( τ ) , . . . , r ( τ T ) } ,where the Kuramoto order parameter is defined as r ( t ) = N (cid:12)(cid:12)(cid:12)(cid:80) Nj =1 e iθ j ( t ) (cid:12)(cid:12)(cid:12) . The Kuramoto order parameter takesthe value of [0 ,
1] and r ( t ) = 1 if the oscillators arein complete synchronization. We used the supervisedlearning scheme to compare the performance between themethod using time-variant topological features and thatusing temporal Kuramoto order parameters. For eachtype of coupling configuration, we randomly split 100 (b)(a) (c) ∞ ∞ ∞ FIG. 3. Examples of time-variant topological features corresponding to the coupling oscillator dynamics in Fig. 2 (a), (b),and (c), respectively. The top row illustrates the three-dimensional persistence diagrams for the loop patterns, while thebottom row illustrates the distribution of the death radii of the connected components appearing along the phase evolutionof the oscillators. As the oscillators approached forming one synchronized-state cluster, the loops quickly disappeared, andthe connected components quickly merged to form one component [Fig. 3(a)]. As the oscillators tended to be divided intomulticlusters of synchronized state, the birth radii of the loops increased, and more connected components survived for a longerperiod of evolution time [Fig. 3(b)(c)]. The number of connected components surviving for a longer period of time correspondedwith the number of clusters that were in the synchronized state. (a) (b)
Time step (T) A cc u r ac y T = T = T = T = Time-variant featuresOrder parameter
FIG. 4. (a) The kernel principal component projection of time-variant topological features obtained from τ = τ , τ , . . . , τ T ,where τ = 0 , τ k +1 − τ k = ∆ τ = 0 . k = 1 , , . . . , T − T is the number of time steps. We illustrate the projection with T = 4 , , ,
11 from left to right. The different colors represent the realization of data with different synchronization schemesas single-cluster (purple), two-cluster (blue), and four-cluster (orange) synchronization. (b) The classification of synchronizedstates using time-variant topological features (blue line) or temporal Kuramoto order parameters (orange line). The linesdepict the average test accuracy across 100 random train-test splits at each value of T . For each split, we randomly split 100realizations into 75 realizations for training and 25 realizations for testing. The shaded areas indicate the confidence intervalsof one standard deviation calculated using the same ensemble of runs. realizations into 75 realizations for training and 25 real-izations for testing. We then applied the support vectormachine [27] for the classification where the input datawere the kernel Gram matrices of the persistence dia-grams or the vectors of the Kuramoto order parameters.Here, we classified the data into three labels that corre- sponded with single-cluster, two-cluster, and four-clustersynchronization, then recorded the accuracy in identify-ing the true labels of the test data.Figure 4(b) depicts the average test accuracy over 100train-test random splits at each value of T . Here, weplotted the topological method using time-variant topo-logical features (blue line), and the order method usingorder parameters (orange line). The shaded area indi-cates the confidence intervals of one standard deviationcalculated using the same ensemble of runs. In general,increasing T served to increase the accuracy level forboth methods since more information on the phase evo-lution was gathered. As Fig. 4(b) shows, the topologicalmethod demonstrated reasonably high accuracy in iden-tifying the synchronized states, even in the early stage ofthe dynamics, i.e., over 85% when T ≥
7. It is clear thatthe corresponding accuracy of the order method was, atmost, 65% at T = 7. As T was increased, the oscilla-tor phases converged into synchronized states, and witha sufficiently large T ( T ≥ B. Chimera states
We applied the time-variant topological features to in-vestigate the evolution of both the coherent and the in-coherent dynamics in oscillatory systems, including thechimera states. A chimera state refers to a state where os-cillators emerge to form two regions of mutually coherentand incoherent populations [28, 29]. Chimera states wereinitially explored in terms of homogeneous oscillator sys-tems [30, 31], before they were then demonstrated exper-imentally [32, 33] to establish their connection with real-world systems such as the human brain networks [34].With reference to Refs. [35, 36], we generated chimerastates for the Kuramoto model by setting the couplingstrength g ij described in Eq. (1) as follows: g ij = π γ , if cos (cid:18) π ( i − j ) N (cid:19) > cos( πη )0 , otherwise . (5)Here, γ is a tunable parameter characterizing the cou-pling strength among the oscillators, and η ∈ [0 ,
1] isa parameter used to control the range of the non-localcoupling. In our numerical simulation, we set N = 256, ω i = 0, α = 1 . η = 0 .
6, and γ = 0 . γ = 6 (for asynchronized states).The coherence-incoherence transition is illustrated bythe time trace of the local order parameter [37, 38] de-fined at each oscillator as follows: l j ( t ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) δ + 1 (cid:88) | j − k |≤ δ e iθ k ( t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (6)where j = 1 , , · · · , N , and θ k represents the phase ofoscillator k in a region of side length 2 δ + 1 centered atoscillator j . The local order parameter quantifies the de-gree of the coherent and incoherent regions around each oscillator and yielded the local properties of the chimerastates. More specifically, oscillator j at time t belongs tothe coherent domain where l j ( t ) ≈
1, or the incoherentdomain where l j ( t ) has much lower values. Figure 5(a)shows three time profile examples for the local order pa-rameters with δ = 12. In the first example, l j ( t ) wereclose to 1 for all j and t >
10, which means it can beconcluded that the oscillators entered global synchroniza-tion as t >
10. In the second example, the values of l j ( t )for 80 ≤ j ≤
200 were lower than 0.6, while the other l j ( t ) were higher than 0 . t >
10. Therefore, thechimera states evolved at t >
10 as the oscillators wereroughly divided into a coherent area and an incoherentarea according to values of l j ( t ). In the final example, theincoherent area dominated since the values of the localparameters were lower than 0.5.The time trace of local order parameters provides auseful indicator for qualitatively evaluating the chimerastates. However, it is not a straightforward task to deter-mine the side length 2 δ +1 of the local region surroundingeach oscillator. The larger the value of δ , the more globalcoherent domains will be captured; however, incoherentdomains will also merge into the global coherent domains.Meanwhile, with a smaller δ value, while more incoherentand spatial domains will be identified, the global coher-ent domains will not be recognized. We thus decided thatthe key aspect of characterizing the chimera states is toexamine the coexistence of two different domains sepa-rated in space, where one part of the oscillator network isoperating coherently while the other exhibits incoherentbehavior. This encouraged us to employ the topologicalmethod to characterize the coherent and incoherent do-mains. By using the time-variant topological features totrack the time trace of these domains, we could betterunderstand how chimera states are evolved, which wouldallow for a qualitative prediction of the chimera states inthe early stage of the dynamics.Here, we present a mapping to transform the j th os-cillator to a point cloud on a torus surface ϕ : θ j → (cid:0) x θ j , y θ j , z θ j (cid:1) , (7)where x θ j = ( R m + R p cos θ j ) cos (cid:18) πjN (cid:19) , (8) y θ j = ( R m + R p cos θ j ) sin (cid:18) πjN (cid:19) , (9) z θ j = R p sin θ j . (10)In our simulations, we set the major radius R m and mi-nor radius R p of the torus as R m = 4 and R p = 1, re-spectively. This mapping demonstrated that we can usehigher-order topological structures such as loops to evalu-ate the chimera states of oscillator systems. For example,in the global synchronized state, the mapped points willtend to be distributed along one major loop on the torussurface, while more minor loops will be formed with theincrease in incoherent regions. (a) (b) (c) (d) FIG. 5. Three examples of different phase dynamics: synchronized state (top row), chimera state (middle row), and asyn-chronized state (bottom row) along the evolution timeline. Their differences can be evaluated by (a) the time trace of localorder parameters, or time-variant topological features such as (b) the three-dimensional persistence diagrams of loop patterns,and (c) the distribution of the death radii of the connected components appearing along the evolution timeline. Persistencediagrams were obtained from the mapped points of the oscillator phases on a torus surface. (d) The shape of the mappedpoints corresponds with each phase of dynamics at t = 40. (a)(b) Synchronized Chimera Asynchronized
FIG. 6. The dimensional reduction of topological similarity features for (a) connected components (top row), and (b) loops(bottom row) via kernel principal component analysis. These figures represent the distribution of each dynamical case when T = 5 , , , , ,
30 from left to right. Each point represents a synchronized (purple), chimera (blue), and asynchronizedstate (orange).
Figures 5(b),(c) present the three-dimensional persis-tence diagrams for the loop patterns and the distribu-tion of the death radii of the connected components ap-pearing along the phase evolution of the oscillators thatcorrespond with the examples in Fig. 5(a). The coloredbar for each plot in Figs. 5(b),(c) indicates the densityof the points in the corresponding plot. The patternsof the mapped points in three-dimensional space corre-sponding with these examples at t = 40 are illustratedin Fig. 5(d). If the oscillators were in global synchro-nized states, for example, the top row of Fig. 5(a) at t ≥
25, only one loop along the torus tube appeared inthe mapped space, meaning we obtained one point in thethree-dimensional diagram in Fig. 5(b) for each t . Asboth coherent and incoherent dynamics emerged in thechimera state in the middle row of Fig. 5, small loopsappeared around the minor circles of the torus surface,which means more points were generated along the timeaxis in the three-dimensional diagram of the loop pat-terns. As illustrated in the bottom row of Fig. 5, the den-sity of these points in the persistence diagram increasedas the incoherent dynamics dominated in the phase dy-namics. The differences among the oscillator phase dy-namics could also be evaluated by observing the distri-bution of the death radii of the connected componentsillustrated in Fig. 5(c). In the global synchronized state,connected components emerged and then quickly mergedat almost the same death radii (top of Fig. 5(c)). Con-versely, in the asynchronized state, the mapped points ofthe oscillators on the torus surface were randomly dis-tributed. The death radii were concentrated in the rangeof 0 to 1, and their distributions were almost the samethroughout time-evolving (bottom of Fig. 5(c)). In thechimera state, the death radii of the connected compo-nents in the coherent region were smaller than those inthe incoherent region, with the death radius exhibitinga wider distribution along the timeline of the chimerastate (middle of Fig. 5(b)). It should be noted thattime-variant topological features provide a novel meansof recognizing chimera states in quantitative terms with-out having to rely on the tuning parameter δ of the localorder parameter.Next, we demonstrate that the kernel method basedon time-variant topological features can also be used tocharacterize the chimera states in early-stage dynamics.Specifically, we prepared 150 cases of temporal phasedata that were differentiated in terms of the initial phasecondition. Here, we considered three labels of synchro-nized, chimera, and asynchronized states, with 50 casesfor each state. We relied on the global Kuramoto or-der parameter r ( t ) at a sufficiently large t to label thedynamics as synchronized states for r ( t = 40) > . . ≤ r ( t = 40) ≤ . r ( t = 40) < .
45. Figure 6 high-lights the projection up to the third component of the kernel principal component analysis for the Gram ma-trix of a total of 150 three-dimensional persistence dia-grams. The purple, blue, and orange points representthe synchronized, chimera, and asynchronized states, re-spectively. To compute the three-dimensional persistencediagrams, we set the time step as τ = 0 , τ k +1 − τ k = 1( k = 1 , , . . . , T − T = 5, it proved difficult to observe any clear dif-ferences among the points presenting different behaviors.As T was increased, the separability increased for bothkernels using the connected components and loops, evenin early stages such as T = 20 , IV. CONCLUDING REMARKS ANDDISCUSSION
In this paper, we demonstrated that the time-varianttopological features constructed from the phase evolu-tion in oscillatory systems can be used to characterizethe behavior of the dynamics, even in the early stagesof the evolution. Such behaviors include global synchro-nization, multicluster synchronization, and chimera stateemergence, which conventional order parameters fail tosufficiently recognize. This indicates that our topologicalapproach is an effective approach for understanding thephase dynamics of oscillators.In previous applications of persistent homology in re-lation to oscillatory systems, only the average temporalpatterns were considered [39, 40]. Our approach fun-damentally differs from such an approach in that it al-lows us to trace the temporal patterns, which are morehelpful to investigating the specific behavior of dynamics.Furthermore, by combining our approach with the ma-chine learning kernel method, we provided an unsuper-vised scheme to characterize the phase dynamics withoutpredefined label training. This aspect is highly signifi-cant from the physical perspective, since unknown dy-namics can be revealed using this unsupervised scheme,including in terms of characterizing the different types ofchimera state.It remains unclear as to whether mapping from a set ofoscillator phases to a point cloud can be regarded as op-timal mapping. In fact, it can be argued that other map-ping methods involving various manifolds could extractmore meaningful and higher dimensional topological in-formation. Moreover, in addition to the values of thephases, other information, such as the phase derivatives,could be used to construct the time-variant topologicalfeatures. In view of this, we expect that our study willbe successfully applied to more practical situations in thefuture, including research involving noisy environments,nonuniform coupling strengths, or asymmetrical networkstructures, all of which may have oscillator networks withtopological configurations that change over time.0 [1] T. Chouzouris, I. Omelchenko, A. Zakharova, J. Hlinka,P. Jiruska, and E. Sch¨oll, Chaos , 045112 (2018).[2] C. J. Stam, Nat. Rev. Neurosci. , 683 (2014).[3] A. Fornito, A. Zalesky, and M. Breakspear, Nat. Rev.Neurosci. , 159 (2015).[4] R. Ferri, S. Miano, O. Bruni, J. Vankova, S. Nevsimalova,S. Vandi, P. Montagna, L. Ferini-Strambi, and G. Plazzi,Clin. Neurophysiol. , 2675 (2005).[5] C. Dibner, U. Schibler, and U. Albrecht, Annu. Rev.Physiol. , 517 (2010).[6] B. A. M. Owens, M. T. Stahl, N. J. Corron, J. N. Blakely,and L. Illing, Chaos , 033109 (2013).[7] N. J. Corron and J. N. Blakely, Proc. Royal Soc. A ,20150222 (2015).[8] S. H. Strogatz, Physica D , 1 (2000).[9] H. Edelsbrunner, D. Letscher, and A. Zomorodian, in Proceedings 41st Annual Symposium on Foundations ofComputer Science (2000) pp. 454–463.[10] Q. H. Tran and Y. Hasegawa, Phys. Rev. E , 032209(2019).[11] A. Myers, E. Munch, and F. A. Khasawneh, Phys. Rev.E , 022314 (2019).[12] H. Kim, J. Hahm, H. Lee, E. Kang, H. Kang, and D. S.Lee, Brain Connect. , 245 (2015).[13] M. Hajij, B. Wang, C. Scheidegger, and P. Rosen, in (2018) pp. 125–134.[14] Q. H. Tran, V. T. Vo, and Y. Hasegawa, Phys. Rev. E , 032308 (2019).[15] Q. H. Tran, M. Chen, and Y. Hasegawa, Preprint atarXiv:2004.03169 (2020).[16] H. Edelsbrunner and J. Harer, Computational Topology.An Introduction (American Mathematical Society, Prov-idence, RI, 2010).[17] J. Reininghaus, S. Huber, U. Bauer, and R. Kwitt, in
Proceedings of the 28th IEEE Conference on ComputerVision and Pattern Recognition (IEEE, Boston, MA,USA, 2015) pp. 4741–4748.[18] G. Kusano, Y. Hiraoka, and K. Fukumizu, in
Proceedingsof The 33rd International Conference on Machine Learn-ing , Proceedings of Machine Learning Research, Vol. 48,edited by M. F. Balcan and K. Q. Weinberger (PMLR,New York, New York, USA, 2016) pp. 2004–2013.[19] M. Carri`ere, M. Cuturi, and S. Oudot, in
Proceedingsof the 34th International Conference on Machine Learn-ing , Proceedings of Machine Learning Research, Vol. 70,edited by D. Precup and Y. W. Teh (PMLR, Interna- tional Convention Centre, Sydney, Australia, 2017) pp.664–673.[20] T. Le and M. Yamada, in
Proceedings of the 32nd In-ternational Conference on Neural Information ProcessingSystems , NIPS’18 (Curran Associates Inc., USA, 2018)pp. 10028–10039.[21] S. Jalan and R. E. Amritkar, Phys. Rev. Lett. , 014101(2003).[22] S. Jalan, R. E. Amritkar, and C.-K. Hu, Phys. Rev. E , 016211 (2005).[23] R. E. Amritkar, S. Jalan, and C.-K. Hu, Phys. Rev. E , 016212 (2005).[24] A. Arenas, A. D´ıaz-Guilera, and C. J. P´erez-Vicente,Phys. Rev. Lett. , 114102 (2006).[25] C. Zhou, L. Zemanov´a, G. Zamora, C. C. Hilgetag, andJ. Kurths, Phys. Rev. Lett. , 238103 (2006).[26] J. Shawe-Taylor, C. K. I. Williams, N. Cristianini, andJ. Kandola, IEEE Trans. Inf. Theory , 2510 (2005).[27] C. M. Bishop, Pattern Recognition and Machine Learn-ing (Springer-Verlag, Berlin, Heidelberg, 2006).[28] O. E. Omel’chenko, Y. L. Maistrenko, and P. A. Tass,Phys. Rev. Lett. , 044105 (2008).[29] M. J. Panaggio and D. M. Abrams, Nonlinearity , R67(2015).[30] Y. Kuramoto and D. Battogtokh, Nonlinear Phenom.Complex Syst. , 380 (2002).[31] D. M. Abrams and S. H. Strogatz, Phys. Rev. Lett. ,174102 (2004).[32] M. R. Tinsley, S. Nkomo, and K. Showalter, Nat. Phys. , 662 (2012).[33] A. M. Hagerstrom, T. E. Murphy, R. Roy, P. H¨ovel,I. Omelchenko, and E. Sch¨oll, Nat. Phys. , 658 (2012).[34] K. Bansal, J. O. Garcia, S. H. Tompson, T. Verstynen,J. M. Vettel, and S. F. Muldoon, Sci. Adv. (2019).[35] Y. Zhu, Z. Zheng, and J. Yang, Phys. Rev. E , 022914(2014).[36] N. Yao, Z.-G. Huang, C. Grebogi, and Y.-C. Lai, Sci.Rep. , 12988 (2015).[37] M. Wolfrum, O. E. Omelchenko, S. Yanchuk, and Y. L.Maistrenko, Chaos , 013112 (2011).[38] I. Omelchenko, Y. Maistrenko, P. H¨ovel, and E. Sch¨oll,Phys. Rev. Lett. , 234102 (2011).[39] B. J. Stolz, H. A. Harrington, and M. A. Porter, Chaos , 047410 (2017).[40] S. Maleti´c, Y. Zhao, and M. Rajkovi´c, Chaos26