Functional Control of Oscillator Networks
Tommaso Menara, Giacomo Baggio, Danielle S. Bassett, Fabio Pasqualetti
FFunctional Control of Oscillator Networks
Tommaso MenaraDepartment of Mechanical Engineering, University of California, Riverside, Riverside, CA92521, USAGiacomo BaggioDepartment of Information Engineering, University of Padova, Padova, 35131, ItalyDanielle S. BassettDepartment of Physics & Astronomy, Department of Bioengineering, Department of Electrical& Systems Engineering, Department of Neurology, Department of Psychiatry, University ofPennsylvania, Philadelphia, PA 19104, USAThe Santa Fe Institute, Santa Fe, NM 87506, USAFabio PasqualettiDepartment of Mechanical Engineering, University of California, Riverside, Riverside, CA92521, [email protected] a r X i v : . [ n li n . AO ] F e b bstract Oscillatory activity is ubiquitous in natural and engineered network systems. Yet, understanding the structure-function relationship in oscillator networks remains an unanswered fundamental question of modern science. Inthis work, we present a method to prescribe exact and robust patterns of functional relations from local networkinteractions of the oscillators. To quantify the behavioral synchrony between agents we introduce the notion of functional patterns , which encode the pairwise relationships between the oscillators’ phases – akin to the Pearsoncorrelation coefficient between two time series. The main contribution of this work is the development of amethod to enforce a desired functional pattern by optimally tailoring the oscillators’ parameters. Importantly,our method is agnostic to the scale at which a network system is studied, computationally efficient, theoreticallysound, and presents an interpretable mapping between the structural principles and the functional implicationsof oscillator networks. As a proof of concept, we apply the proposed method to replicate empirically recordedfunctional relationships from cortical oscillations in a human brain, and to redistribute the active power flow inan electrical grid. Our theory finds applications in the design, analysis, and control of oscillator network systems.
Introduction
Complex coordinated behavior of oscillatory agents determines the purpose, activity, and operation of manynetwork systems in both the natural and artificial domains [11, 14, 44]. Such distributed coordinated behavioracross interconnected components is often encoded in patterns of synchrony [25, 54, 55]. For instance, differentpatterns determine the coordinated motion of orbiting particle systems [50, 65], promote successful mating inpopulations of fireflies [9], regulate the active power flow in electrical grids [18], predict global climate changephenomena [60], and enable numerous cognitive functions in the brain [10, 21]. This rich repertoire of patternsemerges from the properties of the underlying interaction network [56]. Yet, despite its critical role, a completecharacterization of the relationship between collective behaviors and the network’s structural parameters hasremained elusive. This work introduces a direct and interpretable mapping between the structure and thefunction of oscillator networks, and presents a method to enforce the robust emergence of desired functionalrelationships between interconnected components.Controlling the collective motion of interdependent units holds a tremendous potential impact in many real-world applications [58]. To achieve this objective, we abstract rhythmic activity of a system as the outputof a network of diffusively coupled oscillators [3, 17]. While synchronization phenomena in oscillator networkshave been studied extensively (e.g., see Refs. [7, 38, 43, 57]), the development of control methods to imposedesired synchrony across the network units has only recently attracted the attention of the research community[6,20,22,35,42]. Taken together, existing research highlights the importance of controlling distinct configurationsof synchrony, but remains mainly focused on the control of macroscopic synchrony observables.Here, we elevate the ambitions put forth by previous work and provide a simple and elegant framework thatallows to optimally control the spatial organization of the network components and their oscillation frequenciesto achieve desired patterns of synchrony – thus, specifying the oscillators’ pairwise synchrony levels. To accountfor the diversity of interaction types, we study both the case where oscillators can be positively and negativelycoupled, and the case where only positive interactions exist. Finally, we apply our methods to an empiricallyreconstructed neuronal network and a power grid.To quantify the pairwise functional relations between oscillatory units, and inspired by the work in Ref. [4],we define a local correlation metric. Given a pair of phase oscillators i and j with phase trajectories θ i ( t ) and θ j ( t ), respectively, we define the correlation coefficient ρ ij = < cos( θ j ( t ) − θ i ( t )) > t , (1)2here < · > t denotes the average over time. A functional pattern is formally defined as the symmetric matrix R = [ ρ ij ] whose i, j -th entry equals ρ ij . Importantly, a functional pattern explicitly encodes all the localinteractions between oscillators, which would not be possible when considering a global observable (e.g., the orderparameter [3, 59]). It is easy to see that, if two oscillators i and j synchronize after a certain initial transient, ρ ij → i and j become phase-locked (i.e., their phase difference remainsconstant over time), their correlation coefficient converges to some constant value after an initial transient. Ifthe phases of two oscillators i and j evolve independently, their correlation value will remain small over time,i.e., | ρ ij | (cid:28) Analytical Results
The objective of this work is to control the network parameters so that interconnected oscillators exhibit desiredcorrelation values. Due to their rich dynamical repertoire and their wide adoption, we utilize Kuramoto oscillators[1,31]. Specifically, we consider a network G of n oscillators where the i -th oscillator is governed by the dynamicalequation ˙ θ i = ω i + n (cid:88) j =1 A ij sin( θ j − θ i ) , (2)with ω i being its natural frequency and A ij the interconnection weight (i.e., the coupling strength) between oscil-lators i and j . We assume that the network is undirected (i.e., A ij = A ji ) and connected, but the interconnectionscheme can be arbitrary. Further, the connections between oscillators can be either cooperative (i.e., A ij >
0) orcompetitive (i.e., A ij <
0) [27].We focus on phase-locked trajectories of Eq. (2) to achieve desired functional patterns. That is, we aim toselect the weights of the network and the natural frequencies of the oscillators to enforce constant and equalfrequencies ˙ θ i for all i = 1 , . . . , n , so that their phase differences are constant over time. Working towards thisgoal, a few questions arise naturally: Are all functional patterns achievable? Which network configurations allowfor the emergence of desired functional patterns? And, if a certain functional pattern can be achieved, is it robustto perturbations? Surprisingly, we reveal that controlling functional patterns reduces to a convex optimizationproblem, whose solution can be characterized explicitly. We summarize our approach in Fig. 1. Optimal Designs and Interventions as Convex Optimization Problems
To formalize the proposed control problem, let us define the phase difference x ij = θ j − θ i , and let x be the vectorstacking all phase differences with i < j . Rewriting Eq. (2) in matrix form leads to (Materials and Methods):˙ θ = ω − BD ( x ) δ, (3)where B is the (oriented) incidence matrix of the network G , D ( x ) is a diagonal matrix of the sine functions inEq. (2), and δ is a vector collecting all the weights A ij with i < j .Notice that phase-locked solutions to Eq. (3) satisfy ˙ θ ∈ Im( ), where Im( ) denotes the subspace spannedby the vector of all ones. Equivalently, we describe a phase-locked trajectory with an explicit vector of phasedifferences ¯ x . Thus, given desired phase-locked differences ¯ x , we seek for the network parameters that satisfy thefollowing algebraic equation: BD (¯ x ) δ = ω − ω mean (cid:44) ¯ ω, (4)where ω mean is the mean natural frequency of the oscillators (see Materials and Methods). The control problem3 ONTROL original network corrected network parameters
A, ω parameters ˜ A, ˜ ω phases θ phases θ abnormal/undesiredfunctional pattern R desired functionalpattern R oscillators o s c ill a t o r s − c o rr e l a t i o n π π t π π t oscillators o s c ill a t o r s − c o rr e l a t i o n Fig 1.
Network control to enforce a desired functional pattern from an abnormal or undesired one. The leftpanel contains a network of n = 7 oscillators (top left panel, line thickness is proportional to the couplingstrength), whose vector of natural frequencies ¯ ω has zero mean. The phase differences with respect to θ (i.e., θ i − θ ) converge to (cid:2) π π π π π π (cid:3) , as also illustrated in the phases’ evolution from random initial conditions(center left panel, color coded). The bottom left panel depicts the functional pattern R corresponding to suchphase differences over time. The right panel illustrates the same oscillator network after a selection of couplingstrengths and natural frequencies have been tuned (in red, the structural parameters A and ω are adjusted to ˜ A and ˜ ω ) to obtain the phase differences (cid:2) π π π π π π (cid:3) , which encode the desired functional pattern in thebottom right. In this example, we have found the closest set (in the (cid:96) -norm sense) of coupling strengths andnatural frequencies to the original ones that enforce the emergence of the desired pattern. Importantly, only asubset of the original parameters has been modified.that yields a desired functional pattern becomes a convex optimization:min α,β (cid:107) α (cid:107) (cid:63) + (cid:107) β (cid:107) (cid:63) (5)subject to BD (¯ x )( δ + α ) = ¯ ω, ω + β = ¯ ω and ( δ + α ) > α and β are the controllable modifications of the network weights and natural frequencies, respectively,and (cid:107) · (cid:107) (cid:63) denotes any vector norm that one may be interested in minimizing.When feasible, the control problem in Eq. (5) can be solved efficiently due to its convex nature, even forextremely large network systems. However, the physical constraints of specific network instances may restrictaccess to either the network weights or the natural frequencies. Below, we provide a comprehensive investigationof the solutions to Eq. (5) and its variations when the aforementioned constraints arise.4 etworks Without Structural Constraints Can Exhibit Any Feasible FunctionalPattern A functional pattern is a particular square matrix defined by the pairwise correlation coefficients among theoscillators. These coefficients depend only on the phase differences x , which can be uniquely determined bya subset of n − x can be expressed as x = θ − θ = ( θ − θ ) + ( θ − θ ) = x + x . In general, there always exists a( n − n ) / × ( n −
1) matrix P such that x = P x min , where x min is a smallest set of phase differences that canbe used to describe all the remaining ones (see Materials and Methods). Henceforth, we use ¯ x min to indicate adesired phase-locked trajectory, and denote a feasible functional pattern as the matrix R whose entries satisfy ρ ij = cos(¯ x ij ) = cos( f ij (¯ x min )), where ¯ x ij = f ij (¯ x min ) denotes that ¯ x ij is a linear combination of the elements in¯ x min .For a specific ¯ x min , the problem of finding a combination of network weights and natural frequencies thatprescribes the desired functional pattern reads, in its most elementary formulation, asfind δ, ω (6)subject to BD (¯ x min ) δ = ¯ ω. Without any additional constraints, the above problem is always feasible. For instance, a solution can bedetermined by arbitrarily assigning the network weights in δ and then computing the natural frequencies thatsolve Eq. (4). Concurrent Designation of Multiple Functional Patterns
The algebraic core of the control problem in Eq. (6) allows to jointly specify multiple desired phase-lockedsolutions ¯ x ( i )min , i ≥
1, to Eq. (3). Since the constraint in Eq. (6) is a linear system of equations, it can be easilymodified to include multiple desired ¯ x ( i )min as follows: BD (¯ x (1)min ) BD (¯ x (2)min )... δ = ¯ ω ¯ ω ... , which admits a solution δ whenever the known vector [¯ ω ¯ ω . . . ] T on the right-hand side of the equation belongs tothe image of the matrix on the left-hand side. Notice that the above system of equations is obtained by stackingthe n equations BD (¯ x ( i )min ) δ = ¯ ω for each of the additional desired equilibria. We refer the interested reader tothe Supplementary Information and Fig. S1 for a detailed application of this method.Being able to concurrently assign multiple equilibrium patterns is crucial to the investigation and design ofmemory systems [53] – where motifs of activity correspond to distinct memories – and complements previouswork on the search of equilibria in oscillator networks [34]. Condition for the Existence of Functional Patterns in Positive Networks
In practice, not all network systems admit both cooperative and competitive interactions. For instance, negativeinteractions are neither allowed nor physically meaningful in networks of excitatory neurons, in power distributionnetworks (where edge weight denotes conductance and susceptance of a transmission line), and in networks ofmobile robots (where interconnections denote the communication signal strength). Additionally, in some cases,we may not be able to access all the structural parameters. We next address the more challenging problem ofprescribing desired functional patterns in positive networks, where only positive couplings are allowed.5hile the problem of finding natural frequencies to achieve a desired synchronization pattern without alteringthe weights can be addressed by solving Eq. (4), practical situations typically require only the modification of theinterconnection scheme. That is, a local adjustment of the coupling strengths. Formally, we study the problemof controlling the network weights δ > δ (7)subject to BD (¯ x min ) δ = ¯ ω and δ > . The above problem differs from Eq. (6) in that we seek for a combination of positive network weights δ thatsatisfies Eq. (4) without altering the oscillators’ natural frequencies.We first show that the existence of a positive δ that satisfies Eq. (4) is only affected by the sign of the entriesin the diagonal matrix D (¯ x min ), and not by their magnitude. We let B : ,(cid:96) indicate the column of B associatedwith the (cid:96) -th phase difference x ij , and ¯ B be the matrix defined as¯ B : ,(cid:96) = B : ,(cid:96) if sin( f ij (¯ x min )) > , f ij (¯ x min )) = 0 , − B : ,(cid:96) otherwise . Then, we notice that the magnitude of any | sin( f ij (¯ x min )) | can be assimilated by the corresponding weight in δ to define a new variable ¯ δ , and we are left with the linear matrix equality ¯ B ¯ δ = ¯ ω . The latter equation revealsthat the magnitude of the nonzero diagonal entries of D (¯ x min ) does not affect the existence of a solution to theproblem in Eq. (7).A condition for the existence of a positive solution to ¯ B ¯ δ = ¯ ω can be derived by analyzing the projectionsof the natural frequencies ¯ ω on the columns of ¯ B . Notice that requiring positive network weights is equivalentto requiring that there exists a linear combination of the columns of ¯ B with nonnegative coefficients that yields¯ ω . Then, in light of this observation, a sufficient condition for the existence of such a linear combination is asfollows. There exists δ > such that BD (¯ x min ) δ = ¯ ω if there exists a set S of indices of columns of ¯ B satisfying: (i) ¯ B T : ,i ¯ B : ,j ∈ { , − } for all i, j ∈ S with i (cid:54) = j ; (ii) ¯ ω T ¯ B : , S > ; (iii) ¯ ω ∈ Im( ¯ B : , S ) . Due to the structure of ¯ B , the only possible values for ¯ B T : ,i ¯ B : ,j are { , , − } . That is, the angle γ between anytwo columns ¯ B : ,i and ¯ B : ,j is either γ = π , γ = π , or γ = π for ¯ B T : ,i ¯ B : ,j = 1 ,
0, or −
1, respectively. Therefore,requiring (i) is equivalent to ensure that there is an angle γ ≥ π between ¯ B : ,i and ¯ B : ,j . Condition (ii) requires that¯ ω has positive projections on the columns ¯ B : , S . Condition (iii) ensures that the image of S contains ¯ ω . Together,(i), (ii), and (iii) ensure that ¯ ω is contained within the cone generated by the columns ¯ B : , S (see Fig. 2A for aself-contained example). Based on these observations, one can conclude that whenever ¯ ω T ¯ B <
0, the problem inEq. (7) is not feasible, since there cannot exist any linear combination of columns of ¯ B with positive coefficientsthat satisfy the equality constraint. A Graph-theoretic Perspective on Positive Networks
From a graph-theoretic point of view, the condition (i) above prescribes that the columns B : , S correspond tointerconnections that do not share common oscillators or that do not describe directed paths. Recall that adirected Hamiltonian path is a path that visits all the oscillators exactly once. In light of this observation, animmediate corollary to the conditions above is stated next.6 B = − − − ¯ ω = − . . . − − δ δ δ ¯ B : , ¯ B : , ¯ B : , ¯ ω B Hamiltonianpath
Fig 2.
Algebraic and graph-theoretic conditions for the existence of a solution to the problem in Eq. (7). (A)The left side illustrates a simple network of 3 oscillators with adjacency matrix ¯ B and vector of naturalfrequencies ¯ ω . The right side illustrates the cone generated by the columns of ¯ B . In this example, S = { , } satisfies the conditions for the existence of δ > ω is contained within the cone generated by thecolumns ¯ B : , S . (B) The (directed) Hamiltonian path described by the columns of ¯ B : , H , with H = { , } , in thenetwork of panel (A). There exists δ > such that BD (¯ x min ) δ = ¯ ω if (i) the network contains a directed Hamiltonian path; (ii) the columns ¯ B : , H describing such a Hamiltonian path satisfy ¯ ω T ¯ B : , H > . The existence of a Hamiltonian path guarantees that the internal products between columns in ¯ B : , H are all −
1, and the positive projection of ¯ ω on such columns guarantees that δ > Assessing and Enforcing Stability of Functional Patterns
Stability of a functional pattern – that is, small perturbations of initial conditions lead to the same patternasymptotically – is a desirable property in many real-world scenarios. Typically, stability of phase-locked tra-jectories can be assessed by analyzing the linearization of the system. Such linearization is described by theJacobian matrix, whose spectrum determines the stability of the interested trajectories. The Jacobian of theKuramoto oscillators in Eq. (3) can be explicitly written as [17]: J = ∂∂θ ˙ θ = − B diag( { A ij cos( f ij ( x min )) } ) B T (cid:124) (cid:123)(cid:122) (cid:125) L ( x min ) , (8)where L ( x min ) is the Laplacian matrix of the network where the interconnections are weighted by the cosines ofphase differences. That is, the interconnection between nodes i and j reads A ij cos( f ij ( x min )). Thus, the spectrumof the Jacobian matrix evaluated at phase-locked trajectories is equal to the one of the negative Laplacian of acosine-scaled network.The spectrum of the Jacobian in Eq. (8) is known to have one zero eigenvalue (due to rotational symmetry ofthe right-hand side of Eq. (2)) and n − L (¯ x min ) comprises nonpositive off-diagonalentries and strictly positive diagonal entries [17]. Notice that this characterization of the spectrum holds whenever | x ij | < π for all interconnected oscillators i and j . However, in the case of a desired functional pattern includingsome | x ij | > π , stability can be easily enforced when both cooperative and competitive interactions are allowed.By specifying the network weights in δ such that A ij > | x ij | < π and A ij < L is guaranteed to be the Laplacian of a network with positive weights only (see Materials and Methods).Finally, we observe that in the particular case where some | x ij | = π , the cosine-scaled network may become7isconnected, as cos( x ij ) = 0. The Laplacian of a disconnected network has more than one zero eigenvalue,implying that marginal stability of the desired phase differences may occur. Functional Patterns Causing Structural Balance are Unstable
In positive networks, we find that assessing and enforcing stability of a desired functional pattern is less straight-forward. In fact, the matrix L ( x min ) may become a signed Laplacian [64], making the previous analysis inappli-cable. A signed Laplacian contains negative off-diagonal entries when | x ij | < π , and positive off-diagonal entrieswhen | x ij | > π .We exploit the notion of structural balance to derive necessary conditions for the positive definiteness of L ( x min ) and, hence, the stability of the Jacobian. We say that a network is structurally balanced if and onlyif its nodes can be partitioned into 2 sets V and V , such that all competitive interconnections connect nodesin V to nodes in V , and all cooperative interconnections link nodes within V i only to nodes within the samegroup. It has been shown that if a network is structurally balanced, then its Laplacian has mixed eigenvalues [64].Therefore, we can conclude that patterns associated to a structurally balanced cosine-scaled network are unstable.This result highlights the importance of network topology in enabling and determining the stability of certainfunctional patterns. We refer the interested reader to the Supplementary Information for additional results onthe special cases of line and cycle networks. The Supplementary Information also contains a heuristic methodthat posits the lessening or pruning of the interconnections associated with positive off-diagonal entries of theLaplacian to promote stability of functional patterns in positive networks.In summary, the stability of functional patterns can be probed by investigating the spectrum of the Laplacianmatrix of a cosine-scaled network. This network is derived from the desired phase differences and its weightsare defined as A ij cos( f ij (¯ x min )). In positive networks, there is no general condition to test for stability of ¯ x min simply from the analysis of the interconnection topology and network weights. However, if the phase differencessatisfy | x ij | < π , then the Laplacian is well defined and its spectrum is known to guarantee stability of Eq. (8).In the cases where the cosine-scaled network is structurally balanced, the desired pattern is unstable, but thetuning or pruning of the network weights may recover the stability of such a pattern. Applications to Complex Networks
The remainder of this paper contains applications of our theory and methods to an empirically reconstructedbrain network and to the IEEE 39 New England power distribution network. In the former case, we employed theKuramoto model to map structure to function, and found that local metabolic changes underlie the emergenceof functional patterns of recorded neural activity. In the power grid, we utilized our optimization framework torestore the nominal power flow after a fault.
Local Metabolic Changes Govern the Emergence of Distinct Functional Patternsin the Brain
The brain can be studied as a network system in which Kuramoto oscillators approximate the rhythmic activityof different brain regions [10, 15, 28, 35, 45, 46]. Over short time frames, the brain is capable of exhibiting a richrepertoire of functional patterns while the network structure and the interconnection weights remain unaltered.Functional patterns of brain activity not only underlie multiple cognitive processes, but can also be used asbiomarkers for different psychiatric and neurological disorders [49].To shed light on the structure-function relationship of the human brain, we utilize Kuramoto oscillatorsevolving on an empirically reconstructed brain network. We put forth that the intermittent emergence of diversepatterns stems from changes in the oscillators’ natural frequencies – which can be thought of as endogenous8hanges in metabolic regional activity or exogenous interventions to modify undesired synchronization patterns.First, we show that phase-locked trajectories of the Kuramoto model in Eq. (2) can be accurately extracted fromnoisy measurements of neural activity and are an accurate approximation of neural activity time series.We employ structural (i.e., interconnections between brain regions) and functional (i.e., time series of recordedneural activity) data from Ref. [46]. Structural connectivity consists of a sparse weighted matrix A whose entriesrepresent the strength of the physical interconnection between two brain regions. Functional data comprise timeseries of neural activity recorded through functional magnetic resonance imaging (fMRI) of healthy subjects asthey rested. Because the phases of the measured activity have been shown to carry most of the informationcontained in the slow oscillations recorded through fMRI time series, we follow the steps in Ref. [46] to obtainsuch phases from the data (see also Supplementary Information). Next, since frequency synchronization isthought to be a crucial enabler of information exchange between different brain regions and homeostasis of braindynamics [39,40,61], we focus on functional patterns that arise from phase-locked trajectories, as compatible withour analysis. For simplicity, we restrict our study to the cingulo-opercular cognitive system, which comprises n = 12 interacting brain regions [47].To verify whether the Kuramoto model in Eq. (2) is a good approximation of frequency-synchronized neuraldynamics, we identify time windows in the fMRI time series where the signals are phase-locked, and extractthe phases by solving the nonconvex phase synchronization problem [8]. Next, we define the n × n matrix F ofPearson correlation coefficients (also known as functional connectivity) between the filtered time series, and thematrix R = [ ρ ij ] (as in Eq. (1)), which represents the functional pattern in the phase domain. Strikingly, we findthat (cid:107) F − R (cid:107) (cid:28) ω = BD (¯ x min ) δ , where ¯ x min are phase differences obtained from the previous step, and integrate the Kuramotomodel in Eq. (2) with random initial conditions close to ¯ x min . Fig. 3 illustrates our approach to model func-tional connectivity and includes the matrix F for a select time window, along with the matrix R obtained fromintegrating the Kuramoto oscillators. It can be seen that the assignment of natural frequencies according tothe extracted phase differences promotes stable synchronization to accurately replicate the empirical functionalconnectivity F .These results corroborate the postulate that structural connections in the brain support the intermittentactivation of specific functional patterns during rest through regional metabolic changes. Furthermore, we showthat the Kuramoto model represents an accurate and interpretable mapping between the brain anatomical orga-nization and the functional patterns of frequency-synchronized neural co-fluctuations. Once a good mapping isinferred, it can be used to define a generative brain model to replicate in silico how the brain efficiently enactslarge-scale integration of information, and to develop personalized intervention schemes for neurological disordersrelated to abnormal synchronization phenomena [16, 36, 62]. Power Flow Control in Power Networks for Fault Recovery and Prevention
Efficient and robust power delivery in electrical grids is crucial for the correct functioning of this critical infras-tructure. Modern, reconfigurable power networks are expected to be resilient to distributed faults and maliciouscyber-physical attacks [13,63] while being able to rapidly adapt to varying load demands. Therefore, there existsa dire need to design control methods to efficiently operate these networks and react to unforeseen disruptiveevents.The Kuramoto model in Eq. (2) has been shown to be particularly relevant in the modeling of large distributionnetworks and microgrids [18, 51]. Additionally, preliminary work on the control of frequency synchronization in9 rain t b r a i n r e g i o n s t b r a i n r e g i o n s brain regions b r a i n r e g i o n s functional connectivity F − c o rr e l a t i o n brain regions b r a i n r e g i o n s structural connectivity A . m ag n i t ud e brain regions b r a i n r e g i o n s functional pattern R − c o rr e l a t i o n frequency tuning correlation < cos( θ j − θ i ) > t anatomyfiltered fMRItime series ¯ x min Kuramoto
Fig 3.
Conceptual pipeline to reproduce empirically recorded functional connectivity in the brain throughtuning of the natural frequencies of Kuramoto oscillators. The anatomical brain organization provides thenetwork backbone over which the oscillators evolve. The filtered fMRI time series provide the relative phasedifferences between co-fluctuating brain regions, and thus define ¯ x min , which is used to calculate the metabolicchange encoded in the oscillators’ natural frequencies. In this figure, we select the 40-second time window from t = 498 seconds to t f = 538 seconds for subject 18 in the second scanning session. We obtain (cid:107) R − F (cid:107) = 0 . ω i = p (cid:96) i is the active power load at node i , and a ij = | v i || v j | Im( Y ij ), with v i denoting the nodal voltage magnitude and Y ij being the admittance matrix.Notice that, when the phase angles θ of the network nodes are phase-locked and a ij is fixed, the active powerflow is given by a ij sin( θ j − θ i ).We posit that solving the problem in Eq. (5) to design a local reconfiguration of the network parameterscan recover the power distribution before a line fault or provide the minimum parameter tuning to steer theload powers to desired values. In practice, control devices such as Flexible Alternating Current TransmissionSystems (FACTs) allow operators and engineers to change the network parameters [32]. We demonstrate theeffectiveness of our approach by recovering a desired power distribution in the IEEE 39 New England power10 fault
13 245 987611102019 182316 24 15 1413 122221 17 28 29 252627 connection torest of Canada/US B − . − . t − . − . t − . − . t θ buses bu s e s c o rr e l a t i o n buses bu s e s c o rr e l a t i o n buses bu s e s c o rr e l a t i o n n o r m a l f a u l t r e c o v e r y t .
01 50 50 .
01 90 90 . fault control − . − . θ Fig 4.
Fault recovery in the IEEE 39 New England power distribution network through minimal and localintervention. (A) New England power distribution network. The generator terminal buses illustrate the type ofgenerator (coal, nuclear, hydroelectric). We simulate a fault by disconnecting the transmission line 25 (betweenloads 13 and 14). (B) The fault causes the voltage phases θ to depart from normal operating conditions, whichcould cause overheating of some transmission lines (due to violation of the nominal thermal constraint limits)or abnormal power delivery. To recover the pre-fault active power flow, we solve the optimization in Eq. (5) byminimizing the 1-norm of the structural parameter modification δ . The network returns to the initial operativeconditions with a localized modification of the neighboring transmission lines’ impedances.distribution network after a fault. During a regime of normal operation, we simulate a fault by disconnecting abranch between two loads and solve Problem Eq. (5) to find the minimum modification of the couplings a ij thatrecovers the nominal power distribution.We first utilize MATPOWER [66] to solve the power flow problem. Then we use the active powers p (cid:96) and voltages v at the buses to obtain the natural frequencies ¯ ω and the adjacency matrix A of the oscillators,respectively, while the voltage phase angles are used as initial conditions θ (0) for the Kuramoto model in Eq. (2).We integrate the Kuramoto dynamics and let the voltage phases θ ( t ) reach a frequency-synchronized steady state,which corresponds to a normal operating condition. The phase differences also represent a functional patternacross the loads. Next, to simulate a line fault, we disconnect one branch. By solving Problem Eq. (5) with ¯ x min corresponding to the pre-fault steady-state voltage phase differences, we compute the smallest variation of theremaining branch parameters (i.e., admittances) so that the original functional pattern can be recovered. Fig. 4Billustrates the effectiveness of our procedure at recovering the nominal pattern of active power flow by means ofa minimal and localized intervention (see also Fig. S4). Discussion
Distinct configurations of synchrony govern the functions of oscillatory network systems. This work presents asimple and mathematically sound mapping between the structural parameters of arbitrary oscillator networksand their components’ functional interactions. We demonstrate that the control of patterns of synchrony can becast as optimal (convex) design and tuning problems. We also investigate the feasibility of such optimizationsin the cases of networks that admit negative coupling weights and networks that are constrained to positivecouplings. Convexity guarantees that, when feasible, a solution can be found efficiently, even for extremely largenetworks. As a corollary of our control framework, we provide a surprisingly simple method to impose multipledesired equilibria in the phase difference dynamics of Kuramoto oscillators.We remark that our results are also intimately related to the long standing economic problem of enhancing11etwork operations while optimizing wiring costs. In any complex system where synchrony between compo-nents ensures appropriate functions, it is beneficial to maximize synchronization while minimizing the physicalvariations of the interconnection weights [37]. Compatible with this principle, neural systems are thought tohave evolved to maximize information processing by promoting synchronization through optimal spatial organi-zation [2, 12]. Inspired by natural efficiency, Eq. (5) could be utilized to design optimal interaction schemes forlarge-scale computer networks whose performance relies on synchronization-based tasks [30].Prescribing patterns in networks of oscillators has been partially achieved by some authors. In these works,groups of oscillators are forced to synchronize, thus implying that the associated diagonal blocks of the functionalpattern R display values close to 1. Seminal work in Ref. [29] developed a nonlinear feedback control to changethe coupling functions in Eq. (2) to engineer clusters of synchronized oscillators, whereas the authors of Ref. [19]propose the formation of clusters through selective addition of edges to the network. More recently, the controlof partially synchronized states with applications to brain networks is studied in Ref. [35] by means of structuralinterventions, and in Ref. [6] via exogenous stimulation. Perhaps the work that is closest to our approach isRef. [20], where the authors tailor interconnection weights and natural frequencies to achieve a specified level ofphase cohesiveness in a network of heterogeneous Kuramoto oscillators. Our work improves considerably uponthis latter study, whose goal is only to upper bound the phase differences, by enabling the prescription of pairwisephase differences and by investigating the stability properties of all desired functional patterns. All in all, theresults presented in this work go beyond the control of macroscopic synchronization observables and completeprevious research by allowing to specify the synchrony level of pairwise interactions. Methodological Considerations
The framework presented in this work has limitations, which can be addressed in follow up studies. First,despite its capabilities in modeling numerous oscillatory network systems, the Kuramoto model cannot capturethe amplitude of the oscillations, making it most suitable for oscillator systems where most of the informationis conveyed by phase interaction, as demonstrated in Ref. [46] for brain activity at rest. Second, the use ofphase-locked trajectories is instrumental to the control and design of functional patterns. Yet, it is not necessary.In fact, restricting the control to phase-locked dynamics does not capture exotic dynamical regimes in whichonly a part of the network remains frequency-synchronized [33]. A third limitation may arise in situations wherenetwork parameters are not fully known. While being still an active area of research, model reconstruction ofoscillator systems is possible and may be used in such scenarios [41]. Fourth and finally, albeit not explicitlyconsidered in this work, sparsity constraints can also be added to Eq. (5) to study cases where only a subset ofthe network weights is accessible for modification (i.e., δ consists only of a subset of the nonzero entries in A ). Conclusion
In this work, we have provided a method to achieve desired functional patterns in complex oscillator networks.Specifically, we have demonstrated that the structural network parameters that yield desired patterns of syn-chrony of the oscillators’ phases can be obtained as the solution to convex optimization problems, and that thestability of such patterns can be probed by investigating the spectrum of a Laplacian matrix. To showcase theapplicability of our framework, we have replicated functional patterns that arise from brain activity at rest, andrecovered the nominal active power flow pattern in an electrical grid after a simulated fault. Our results shedlight on the relationship between network structure and function in all domains where synchronization plays acrucial role, such as the human brain and numerous engineered network systems.Directions of future research can be both of a theoretical and practical nature. For instance, follow up studiescan focus on the derivation of a general condition to check for the stability of any feasible functional pattern inpositive networks. Further, a thorough investigation of which network structures allow for the prescription of12ultiple equilibria may be particularly relevant in the context of memory systems, where different patterns areassociated to different memory states. Specific practical applications may also require the inclusion of sparsityconstraints on the accessible structural parameters for the implementation of the proposed control and designframework.
Materials and Methods
Matrix Form of Eq. Eq. (2) and Phase-locked Solutions
For a given network of oscillators, we let the entries of the (oriented) incidence matrix B be defined component-wise after choosing the orientation of each interconnection ( i, j ) such that i points to j if i < j : B k(cid:96) = − k is the source node of the interconnection (cid:96) , B k(cid:96) = 1 if oscillator k is the sink node of the interconnection (cid:96) , and B k(cid:96) = 0 otherwise. Then, the matrix form of Eq. Eq. (2) can be written as˙ θ = ω − B . . . sin( x ij ) . . . δ = ω − B diag( { sin( x ij ) } ) δ == ω − BD ( x ) δ, where D ( x ) is the diagonal matrix of the sine functions in Eq. Eq. (2).Phase-locked solutions to the above equation (i.e., Eq. Eq. (3)) satisfy ˙ θ = k , where k is the mean naturalfrequency. Since T B = 0, the only k such that ( ω − k ) belongs to the image of B is k = ω mean = n (cid:80) ni =1 ω i .Notice that, from an algebraic point of view, the operation ¯ ω = ω − k is equivalent to projecting ω onto theorthogonal subspace to . That is, ¯ ω satisfies T ¯ ω = 0. Any Feasible Functional Pattern has n − Degrees of Freedom
The values of a functional pattern can be uniquely specified using a set of n − x = M θ , where M is a ( n − n ) / × n matrix whose k -th row, associated to x ij , is all zeros except for b ki = − b kj = 1. Consider the first n − M , associated to x , x , . . . x n ,and notice that they are linearly independent. Moreover, the row associated to x ij , i >
1, can be obtained bysubtracting the row associated to x i to the row associated to x j , implying that the row rank of M is n − n − M defines a full row-rank matrix M min (e.g., any n − x min = M min θ , where x min isa smallest set of phase differences that can be used to quantify synchronization among all oscillators. Becauseker( M min ) = , it holds that M † min x min = θ + c , where M † min denotes the Moore-Penrose pseudo-inverse of M min and c is any real number. The latter equality implies that we can obtain the phases θ from x min modulo rotation,and yields M M † min x min = M ( θ + c ) = M θ + 0 = x. The above equation reveals that all the differences x are encoded in x min . That is, any x ij can be writtenas a linear combination of elements in x min . For example, if n = 3 and x min = [ x x ] T , then x is a linearcombination of the differences in x min , i.e., x = (cid:104) − − − (cid:105) (cid:2) − − (cid:3) † x min , in which x = x + x . Thus,because n − nforcing Stability of Functional Patterns To ensure that the Jacobian matrix in Eq. Eq. (8) is the Laplacian of a network with positive weights only, wesolve the problem in Eq. Eq. (7) with a slight modification of the constraints. To be specific, we post-multiplythe matrix B in Eq. Eq. (7) as B ∆, where ∆ = diag( { sign(cos( x ij )) } ) is a matrix that changes the sign of thecolumns of B associated to negative weights in the cosine-scaled network. By solving for positive interconnectionweights the problem in Eq. Eq. (7), we enforce a stable Jacobian in a network where the final couplings are ∆ δ . Materials and data availability
Brain data used in the preparation of this work were obtained from the supplementary material of [46] and canbe freely downloaded at https://doi.org/10.1371/journal.pcbi.1004100.s006 . The IEEE 39 New Englanddata parameters and interconnection scheme can be found in the reference textbook [48] (see also SupplementaryInformation). Source code and documentation for the numerical simulations presented here are freely availablein GitHub at: https://github.com/tommasomenara/functional_control . References [1] Juan A Acebr´on, Luis L Bonilla, Conrad J P´erez Vicente, F´elix Ritort, and Renato Spigler. The Kuramotomodel: A simple paradigm for synchronization phenomena.
Reviews of modern physics , 77(1):137, 2005.[2] S. Achard and E. Bullmore. Efficiency and cost of economical brain functional networks.
PLOS Computa-tional Biology , 3(2):e17, 2007.[3] A. Arenas, A. D´ıaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou. Synchronization in complex networks.
Physics Reports , 469(3):93–153, 2008.[4] A. Arenas, A. D´ıaz-Guilera, and C. J. P´erez-Vicente. Synchronization reveals topological scales in complexnetworks.
Physical Review Letters , 96:114102, Mar 2006.[5] T. Athay, R. Podmore, and S. Virmani. A practical method for the direct analysis of transient stability.
IEEE Transactions on Power Apparatus and Systems , 98(2):573–584, 1979.[6] K. Bansal, J. O. Garcia, S. H. Tompson, T. Verstynen, J. M. Vettel, and S. F. Muldoon. Cognitive chimerastates in human brain networks.
Science Advances , 5(4), 2019.[7] I. Belykh, R. Jeter, and V. Belykh. Foot force models of crowd dynamics on a wobbly bridge.
ScienceAdvances , 3(11):e1701512, 2017.[8] N. Boumal. Nonconvex phase synchronization.
SIAM Journal on Optimization , 26(4):2355–2377, 2016.[9] J. Buck and E. Buck. Biology of synchronous flashing of fireflies.
Nature , 211(5049):562–564, 1966.[10] J. Cabral, E. Hugues, O. Sporns, and G. Deco. Role of local network oscillations in resting-state functionalconnectivity.
NeuroImage , 57(1):130–139, 2011.[11] A. Cavagna, A. Cimarelli, I. Giardina, G. Parisi, R. Santagati, F. Stefanini, and M. Viale. Scale-freecorrelations in starling flocks.
Proceedings of the National Academy of Sciences , 107(26):11865–11870, 2010.[12] B. L. Chen, D. H. Hall, and D. B. Chklovskii. Wiring optimization can relate neuronal structure andfunction.
Proceedings of the National Academy of Sciences , 103(12):4723–4728, 2006.[13] S. P. Cornelius, W. L. Kath, and A. E. Motter. Realistic control of network dynamics.
Nature Communi-cations , 4, 2013. 1414] M. C. Cross and P. C. Hohenberg. Pattern formation outside of equilibrium.
Reviews of Modern Physics ,65(3):851, 1993.[15] D. Cumin and C. P. Unsworth. Generalising the Kuramoto model for the study of neuronal synchronisationin the brain.
Physica D: Nonlinear Phenomena , 226(2):181–196, 2007.[16] G. Deco, J. Cruzat, J. Cabral, E. Tagliazucchi, H. Laufs, N. K. Logothetis, and M. L. Kringelbach. Awak-ening: Predicting external stimulation to force transitions between different brain states.
Proceedings of theNational Academy of Sciences , 116(36):18088–18097, 2019.[17] F. D¨orfler and F. Bullo. Synchronization in complex networks of phase oscillators: A survey.
Automatica ,50(6):1539–1564, 2014.[18] F. D¨orfler, M. Chertkov, and F. Bullo. Synchronization in complex oscillator networks and smart grids.
Proceedings of the National Academy of Sciences , 110(6):2005–2010, 2013.[19] R. M. D’Souza and M. Mitzenmacher. Local cluster aggregation models of explosive percolation.
PhysicalReview Letters , 104:195702, May 2010.[20] M. Fazlyab, F. D¨orfler, and V. M. Preciado. Optimal network design for synchronization of coupled oscilla-tors.
Automatica , 84:181–189, 2017.[21] Juergen Fell and Nikolai Axmacher. The role of phase synchronization in memory processes.
Nature ReviewsNeuroscience , 12(2):105–118, 2011.[22] A. Forrow, F. G. Woodhouse, and J. Dunkel. Functional control of network dynamics using designedLaplacian spectra.
Physical Review X , 8(4):041043, 2018.[23] S. Gerˇsgorin. ¨Uber die abgrenzung der eigenwerte einer matrix.
Bulletin de l’Acad´emie des Sciences del’URSS. Classe des sciences math´ematiques et na , pages 749–754, 1931.[24] C. Godsil and G. F. Royle.
Algebraic Graph Theory . Graduate Texts in Mathematics. Springer New York,2001.[25] Martin Golubitsky, Ian Stewart, and Andrei T¨or¨ok. Patterns of synchrony in coupled cell networks withmultiple arrows.
SIAM Journal on Applied Dynamical Systems , 4(1):78–100, 2005.[26] M. Grant, S. Boyd, and Y. Ye. CVX: Matlab software for disciplined convex programming, 2009.[27] H. Hong and S. H. Strogatz. Kuramoto model of coupled oscillators with positive and negative couplingparameters: An example of conformist and contrarian oscillators.
Physical Review Letters , 106(5):054102,2011.[28] F. C Hoppensteadt and E. M Izhikevich.
Weakly Connected Neural Networks . Springer Science & BusinessMedia, 2012.[29] Istv´an Z. Kiss, Craig G. Rusin, Hiroshi Kori, and John L. Hudson. Engineering complex dynamical struc-tures: Sequential patterns and desynchronization.
Science , 316(5833):1886–1889, 2007.[30] G. Korniss, M. A. Novotny, H. Guclu, Z. Toroczkai, and P. A. Rikvold. Suppressing roughness of virtualtimes in parallel discrete-event simulations.
Science , 299(5607):677–679, 2003.[31] Y. Kuramoto. Self-entrainment of a population of coupled non-linear oscillators. In H. Araki, editor,
Int.Symposium on Mathematical Problems in Theoretical Physics , volume 39 of
Lecture Notes in Physics , pages420–422. Springer, 1975. 1532] A. L’Abbate, G. Migliavacca, U. H¨ager, C. Rehtanz, S. R¨uberg, H. Ferreira, G. Fulli, and A. Purvins.The role of FACTS and HVDC in the future paneuropean transmission system development. In , pages 1–8, 2010.[33] Matthew H Matheny, Jeffrey Emenheiser, Warren Fon, Airlie Chapman, Anastasiya Salova, Martin Rohden,Jarvis Li, Mathias Hudoba de Badyn, M´arton P´osfai, Leonardo Duenas-Osorio, et al. Exotic states in asimple network of nanoelectromechanical oscillators.
Science , 363(6431), 2019.[34] D. Mehta, N. S. Daleo, F. D¨orfler, and J. D. Hauenstein. Algebraic geometrization of the Kuramoto model:Equilibria and stability analysis.
Chaos: An Interdisciplinary Journal of Nonlinear Science , 25(5):053103,2015.[35] T. Menara, G. Baggio, D. S. Bassett, and F. Pasqualetti. A framework to control functional connectivity inthe human brain. In
IEEE Conf. on Decision and Control , pages 4697–4704, Nice, France, December 2019.[36] T. Menara, G. Lisi, F. Pasqualetti, and A. Cortese. Brain network dynamics fingerprints are resilient todata heterogeneity.
Journal of Neural Engineering , 2020.[37] T. Nishikawa and A. E. Motter. Maximum performance at minimum cost in network synchronization.
Physica D: Nonlinear Phenomena , 224(1-2):77–89, 2006.[38] T. Nishikawa and A. E. Motter. Network synchronization landscape reveals compensatory structures, quan-tization, and the positive effect of negative interactions.
Proceedings of the National Academy of Sciences ,107(23):10342–10347, 2010.[39] R. Noori, D. Park, J. D. Griffiths, S. Bells, P. W. Frankland, D. Mabbott, and J. Lefebvre. Activity-dependent myelination: A glial mechanism of oscillatory self-organization in large-scale brain networks.
Proceedings of the National Academy of Sciences , 117(24):13227–13237, 2020.[40] Agostina Palmigiano, Theo Geisel, Fred Wolf, and Demian Battaglia. Flexible information routing bytransient synchrony.
Nature Neuroscience , 20(7):1014, 2017.[41] M. J. Panaggio, M. Ciocanel, L. Lazarus, C. M. Topaz, and B. Xu. Model reconstruction from temporal datafor coupled oscillator networks.
Chaos: An Interdisciplinary Journal of Nonlinear Science , 29(10):103116,2019.[42] L. Papadopoulos, J. Z. Kim, J. Kurths, and D. S. Bassett. Development of structural correlations andsynchronization from adaptive rewiring in networks of Kuramoto oscillators.
Chaos: An InterdisciplinaryJournal of Nonlinear Science , 27(7):073115, 2017.[43] Louis M. Pecora and Thomas L. Carroll. Master stability functions for synchronized coupled systems.
Physical Review Letters , 80:2109–2112, Mar 1998.[44] A. Pikovsky, M. Rosenblum, and J. Kurths.
Synchronization: A Universal Concept in Nonlinear Sciences .Cambridge University Press, 2003.[45] A. Politi and M. Rosenblum. Equivalence of phase-oscillator and integrate-and-fire models.
Physical ReviewE , 91(4):042916, 2015.[46] A. Ponce-Alvarez, G. Deco, P. Hagmann, G. L. Romani, D. Mantini, and M. Corbetta. Resting-state tem-poral synchronization networks emerge from connectivity topology and heterogeneity.
PLoS ComputationalBiology , 11(2):1–23, 02 2015. 1647] J. D. Power, A. L. Cohen, S. M. Nelson, G. S. Wig, K. A. Barnes, J. A. Church, A. C. Vogel, T. O. Laumann,F. M. Miezin, B. L. Schlaggar, and S. E. Petersen. Functional network organization of the human brain.
Neuron , 72(4):665–678, 2011.[48] P. W. Sauer and M. A. Pai.
Power System Dynamics and Stability . Prentice Hall, 1998.[49] A. Schnitzler and J. Gross. Normal and pathological oscillatory communication in the brain.
Nature ReviewsNeuroscience , 6(4):285, 2005.[50] R. Sepulchre, D. A. Paley, and N. E. Leonard. Stabilization of planar collective motion: All-to-all commu-nication.
IEEE Transactions on Automatic Control , 52(5):811–824, 2007.[51] J. W. Simpson-Porco, F. D¨orfler, and F. Bullo. Droop-controlled inverters are Kuramoto oscillators.
IFACProceedings Volumes , 45(26):264 – 269, 2012. 3rd IFAC Workshop on Distributed Estimation and Controlin Networked Systems.[52] P. S. Skardal and A. Arenas. Control of coupled oscillator networks with application to microgrid technolo-gies.
Science Advances , 1(7):e1500339, 2015.[53] P. S. Skardal and A. Arenas. Memory selection and information switching in oscillator networks withhigher-order interactions.
Journal of Physics: Complexity , 2(1):015003, dec 2020.[54] F. Sorrentino and E. Ott. Network synchronization of groups.
Physical Review E , 76(5):056114, 2007.[55] I. Stewart, M. Golubitsky, and M. Pivato. Symmetry groupoids and patterns of synchrony in coupled cellnetworks.
SIAM Journal on Applied Dynamical Systems , 2(4):609–646, 2003.[56] S. H. Strogatz. Exploring complex networks.
Nature , 410(6825):268–276, 2001.[57] S. H. Strogatz.
SYNC: The Emerging Science of Spontaneous Order . Hyperion, 2003.[58] S. H. Strogatz and I. Stewart. Coupled oscillators and biological synchronization.
Scientific American ,269(6):102–109, 1993.[59] T. Tanaka and T. Aoyagi. Optimal weighted networks of phase oscillators for synchronization.
PhysicalReview E , 78(4):046210, 2008.[60] K. E. Trenberth. Spatial and temporal variations of the Southern Oscillation.
Quarterly Journal of theRoyal Meteorological Society , 102(433):639–653, 1976.[61] F. Varela, J. P. Lachaux, E. Rodriguez, and J. Martinerie. The brainweb: Phase synchronization andlarge-scale integration.
Nature Reviews Neuroscience , 2(4):229–239, 2001.[62] F. V´aˇsa, M. Shanahan, P. J. Hellyer, G. Scott, J. Cabral, and R. Leech. Effects of lesions on synchrony andmetastability in cortical networks.
Neuroimage , 118:456–467, 2015.[63] J.-W. Wang and L.-L. Rong. Cascade-based attack vulnerability on the US power grid.
Safety Science ,47(10):1332 – 1336, 2009.[64] D. Zelazo and M. B¨urger. On the robustness of uncertain consensus networks.
IEEE Transactions on Controlof Network Systems , 4(2):170–178, 2017.[65] F. Zhang and N. E. Leonard. Coordinated patterns of unit speed particles on a closed curve.
Systems &Control Letters , 56(6):397–407, 2007.[66] R. D. Zimmerman, C. E. Murillo-S´anchez, and R. J. Thomas. MATPOWER: Steady-state operations,planning, and analysis tools for power systems research and education.
IEEE Transactions on power systems ,26(1):12–19, 2010. 17 upplementary Information
Specification of Multiple Phase-locked Equilibria By Tailoring of the NetworkStructural Parameters
The convex optimization problem proposed in the main text allows to tailor the network weights and the naturalfrequencies of the oscillators to specify multiple equilibria for the dynamics of the phase differences x . Theseequilibria correspond to phase trajectories θ that evolve with constant, desired phase differences.We now provide an example where we jointly impose, for a complete graph of n = 7 oscillators, two equilibriafor the phase difference dynamics. Specifically, by taking θ as a reference, we choose two points for the phasedifferences x i = θ i − θ to be set as equilibria: ¯ x (1)min = (cid:2) π π π π π π (cid:3) T and ¯ x (2)min = (cid:2) π π π π π π (cid:3) T . Theinitial network parameters (adjacency matrix and zero-mean natural frequencies) read as: A = and ¯ ω = . − . . . − . − . . , respectively. To impose the desired equilibria, we numerically solve the following convex problem through stan-dard cvx routines [26]: min α (cid:107) δ + α (cid:107) subject to BD (cid:16) ¯ x (1)min (cid:17) BD (cid:16) ¯ x (2)min (cid:17) ( δ + α ) = (cid:34) ¯ ω ¯ ω (cid:35) . The solution δ yields the following corrected adjacency matrix:˜ A = . . − . − . . − . . − . . . . . − . . . . . − . . . − . . . . . − . . . . An investigation of the Jacobian spectrum (see main text for results on stability) of the phase differences computedat the two equilibria ¯ x (1)min and ¯ x (2)min reveals that the first equilibrium point is unstable and that the second oneis locally stable. We illustrate the outcome of the above procedure to specify multiple equilibria in Fig. 5, wherethe phase differences start at ¯ x (1)min at time t = 0, and converge to ¯ x (2)min after a perturbation is applied at t = 50to force them out from the equilibrium ¯ x (1)min . 18 tability Results For Functional Patterns With Angle Differences Larger than π inthe Case of Line and Cycle with Positive Weights Consider a line network of n (ordered) oscillators with positive-only weights that possesses an equilibrium forthe phase difference dynamics satisfying | x i,i +1 | > π for some i ∈ { , . . . , n − } . It is straightforward to deducefrom the results on structural balance [64] that the considered equilibrium is unstable. This implies that thefunctional pattern associated with that equilibrium is unstable.Intuitively, the simplest topology that lends itself to a characterization of stable phase configurations including | x i,j | > π and allowing only positive weights is the cycle (i.e., a line network where the first and last oscillatorsare connected). Consider a cycle network of n > Theorem 1 ( Stability of phase differences equilibria with | x ij | > π in cycle networks ) The equilibrium ¯ x cycle = [ x x . . . x n − ,n ] T = [ γ ϕ . . . ϕ n − ] T with | γ | > π , | ϕ i | < π , and ϕ i = ϕ n − i for all i = 1 , . . . , n − ,is stable if and only if (i) A i,i +1 = A
12 sin( γ )sin( ϕ i − ) for i = 2 , . . . , n , with n − i + 3 (cid:44) if i = 2 , and n = 1 (cid:44) ; (ii) | cot γ | ≤ (tan( ϕ ) + · · · + tan( ϕ n − )) − .Moreover, if ϕ = · · · = ϕ n − and n → ∞ , the largest possible value for γ tends to the value γ ≈ . .Proof of Theorem 1: To assess the stability of the equilibrium ¯ x cycle = [ x x . . . x n − ,n ] T = [ γ ϕ . . . ϕ n − ] T with | γ | > π , | ϕ i | < π , and ϕ i = ϕ n − i for all i = 1 , . . . , n −
1, we analyze the spectrum of the Jacobian J (¯ x cycle ) = −L (¯ x cycle ). From Ref. [64, Corollary IV.7], we have that a necessary and sufficient condition for the Laplacianmatrix L (¯ x cycle ) of the cosine-scaled network to be positive semidefinite is | a cos( γ ) | ≤ R − , (9)with R being the effective resistance of the graph in which the edge (1 ,
2) has been removed. That is, R = 1 a cos( ϕ ) + · · · + 1 a n − ,n cos( ϕ n − ) . (10)Since the adjacency matrix satisfies A = A T and ¯ x cycle is an equilibrium for the difference dynamics of the cyclenetwork, the network weights must be identical pairwise: A i,i +1 = A n − i +2 ,n − i +3 = A sin( γ )sin( ϕ i − ) , (11)for i = 2 , . . . , n with the convention n − i + 3 (cid:44) i = 2. Thus, the angles being identical pairwise implies thatthe network weights must also be identical pairwise, which yields condition (i) of Theorem 1. Moreover, pluggingthe network weights from Eq. (11) into Eq. (12) yields R = tan( ϕ ) a sin( γ ) + · · · + tan( ϕ n − ) a sin( γ ) , (12)which makes the condition in Eq. (9) become, after simple calculations, condition (ii) of Theorem 1: | cot γ | ≤ (tan( ϕ ) + · · · + tan( ϕ n − )) − . (13)19hus, given that the Laplacian L (¯ x cycle ) is positive definite, the Jacobian J (¯ x cycle ) is stable, and its only zeroeigenvalue is due to rotational symmetry of the right-hand side of the Kuramoto dynamics in Eq. (2) of the maintext. The same reasoning can be repeated for n even, and it is therefore omitted here in the interest of space.For ϕ = · · · = ϕ n − , we have that ϕ i = 2 π − γn − . Hence, the right-hand side of Eq. (13) becomes cot( ϕ ) n − , and lim n →∞ cot( ϕ ) n − = π − γ . Since | γ | > π , plugging thelimit value for π − γ into Eq. (13) and solving for the equality yields γ − tan( γ ) = 2 π , whose unique solution is γ ≈ . (cid:4) A Heuristic Method To Promote Stability of Functional Patterns in PositiveNetworks
Although there is no general stability condition for arbitrary positive networks and desired functional patterns, weprovide a heuristic procedure to promote the stability of functional patterns that include negative correlations ina network with nonnegative weights. We recall the definition of Gerschgorin disks and the Gerschgorin Theorem.
Definition of Gerschgorin disk . Let M be a n × n complex matrix. Let the radius of the i -th disk, i = 1 , . . . , n ,be r i = (cid:80) j (cid:54) = i | M ij | . and let the center of the disk be M ii . The disk D i = ( M ii , r i ) is the i -th Gerschgorin disk. Theorem 2 ( Gerschgorin [23]) All the eigenvalues of the matrix M lie within the union (cid:83) ni =1 D i of Gerschgorindisks.To promote stability of functional patterns with phase differences | x ij | > π , it is advantageous to design thenetwork weights by minimizing the ones associated to a cos( x ij ) < A ij as much as possible). Re-ducing the magnitude of these coupling strengths, or even pruning such interconnections, causes the Gerschgorindisks of the Laplacian L ( x min ) to lie almost entirely in the right half-plane. In practice, one may jointly optimizethe network weights to satisfy the equilibrium constraints and the heuristic stability strategy. We let P (resp., N ) denote the set of indices associated with A ij cos(¯ x ij ) > A ij cos(¯ x ij ) < α c (cid:107) α ( P ) (cid:107) (cid:63) + c (cid:107) δ ( N ) + α ( N ) (cid:107) (cid:63) (14)subject to BD (¯ x min )( δ + α ) = ¯ ω, ( δ + α ) > , where (cid:107) · (cid:107) (cid:63) is a desired vector norm, c , c > α ( P ) denotes the entries of the tuning vector α that are associated to positive weights in the cosine-scaled network, and α ( N ) denotes the entries of the tuningvector α that are associated to negative weights δ ( N ) in the cosine-scaled network.To showcase the effectiveness of the proposed heuristic method to promote stability, we construct an examplewith n = 7 oscillators and desired minimum vector of phase differences¯ x (1)min = (cid:20) π π π π π π (cid:21) T , where x ij = θ j − θ , j = 2 , . . . ,
7. Notice that the first difference x > π/
2, hence cos( x ) <
0. Consider the20scillator network with structural parameters that read as: A = . . . . . . . . . . . . and ¯ ω = − . . . − . − . − . . . (15)Such a network admits the following phase-locked equilibrium¯ x (0)min = (cid:104) π π π π π π (cid:105) T , which differs from ¯ x (1)min only in x , and has all x ij satisfying | x ij | < π/ x (0)min and the functional pattern R associated with this equilibrium.The network considered in this example has 9 network weights that can be modified, and the desired equilib-rium ¯ x (1)min implies that N = { } and P = { , , . . . , } . To compute the optimal tuning of the network weights wesolve the optimization in Eq. (14) (through standard cvx routines [26]) by minimizing the (cid:96) -norm. The adjustednetwork adjacency matrix becomes:˜ A = . . . . . . . . . . . . , (16)where the optimal network correction α has disconnected the entry A (highlighted in red). The matrix ˜ A guarantees stability of the functional pattern associated with ¯ x (1)min , as we illustrate in Fig. 6C-D. Phase-locked time series from fMRI data to approximate with Kuramotooscillators
Following the procedure in Ref. [46], we apply a narrow-band filter in the low frequency range [0 .
04 0 .
07] Hzto the time series for each brain region. Next, to obtain a measure of the functional synchrony between thebrain regions, we generate the n × n functional connectivity matrix F , whose entry F ij indicates the pairwisePearson correlation coefficient between filtered time series of recorded neural activity. To map these functionalcorrelations to the phase domain, we extracted the phase time series ˜ θ ( t ) by applying a Hilbert Transform to thefiltered signals. From ˜ θ ( t ), one can identify time windows over which frequency synchronization emerges withthe aid of the phase-locking value matrix P = [ P ij ], where P ij ( t , t f ) = 1 t f − t t f (cid:88) t = t (cid:12)(cid:12)(cid:12) e i [ θ j ( t ) − θ i ( t )] (cid:12)(cid:12)(cid:12) . Clearly, if P ij ( t , t f ) ≈ i, j , then the time window [ t t f ] comprises phase-locked trajectories.Since the phase time series ˜ θ ( t ) are derived from inherently noisy measurements, we compute the best estimate21f the phases θ ∗ (modulo rotation) that are compatible with the noisy measurement in ˜ θ ( t ) by solving the nonconvex phase synchronization problem [8] – that is, the estimation of phases from noisy pairwise relativephase measurements. Given a time window of frequency-synchronized phase time series ˜ θ , we find that R ≈ F (see main text and Fig. 7). Moreover, it holds that (cid:107) R − F (cid:107) → P ij → θ ∗ ( t ) (encoded in the matrix R ) represent the same functionalrelationships that are measured in fMRI data (encoded in the matrix F ), supporting the usage of Kuramotooscillators to analyze remote synchronization. Power network modeling and assumptions
The main manuscript contains an application of our network tuning methods to the IEEE 39 New England powerdistribution network [5, 48]. To model the dynamics of this network, we consider a connected power networkwith generators V and load buses V . A structure-preserving power network model contains |V | second-orderNewtonian and |V | first-order kinematic phase oscillators obeying [48]: M i ¨ θ i + D i ˙ θ i = ω i + (cid:80) |V | j =1 a ij sin( θ j − θ i ) , i ∈ V ,D i ˙ θ i = ω i + (cid:80) |V | j =1 a ij sin( θ j − θ i ) , i ∈ V , (17)where M i , D i are the generator inertia constant, and the damping coefficient, respectively. In the equation forthe generators V , ω i = P m,i , which is the mechanical power input from the prime mover, and in the equationfor the load buses V , ω i = P (cid:96),i , which denotes the real power drawn by load i . Finally, the weight a ij equals a ij = | v i || v j | Im( Y ij ), with v i denoting the nodal voltage magnitude and Y ij being the admittance matrix. Theabove structure-preserving power network model represents an AC grid with a synchronous generator.Owing to Ref. [17, Lemma 1], the existence and local exponential stability of synchronized solutions of theoscillator model Eq. (17) can be entirely described by means of the first-order Kuramoto model. That is, the loaddynamics of a structure-preserving power grid model has the same stable synchronization manifold of Eq. (2) inthe main text.We assume that thermal limit constraints are equivalent to phase cohesiveness requirements. To be precise,we obtain a bounded power flow a ij sin( θ j − θ i ) for the line ( i, j ) whenever the angular distance | θ j − θ i | isbounded, which is satisfied by frequency-synchronized phase trajectories. Moreover, we assume constant voltagemagnitudes | v i | at the loads, so that the weights a ij can be considered fixed. This is a standard assumption inpower systems ( decoupling assumption ). We refer the interested reader to Ref. [17, Remark 1] for further details.The generators and bus parameters for the IEEE New England Power network are available in the originalarticle and in classic textbooks [5, 48]. In our simulations, we utilized the standard optimal power flow solverprovided by MATPOWER to compute the parameters v , p (cid:96) and θ (0) needed to integrate the Kuramoto modelin Matlab through a standard ODE45 solver. 22 erturbation0 25 50 75 100 125 1500 π/ π/ π/ π/ π/ tx x x x x x Fig 5.
Multiple equilibria for the phase difference dynamics. The trajectories start at the unstable equilibriumpoint ¯ x (1) = (cid:2) π π π π π π (cid:3) T at time t = 0, and converge to the locally stable point ¯ x (2) = (cid:2) π π π π π π (cid:3) T after a small perturbation 0 . · p , with p ∈ [0 1] , is applied to the phase difference trajectories at time t = 50. A . . . . tx x x x x x B oscillators o s c ill a t o r s R − m ag n i t ud e C . . . . tx x x x x x D oscillators o s c ill a t o r s R − m ag n i t ud e Fig 6.
Results of the heuristic method promoting stability of functional patterns containing negativecorrelations. (A) Phase differences of the network with adjacency matrix A in Eq. Eq. (15). The phasedifferences converge to ¯ x (1)min . (B) The functional pattern associated with the phase differences in panel (A). (C)Phase differences of the network with adjacency matrix ˜ A in Eq. Eq. (16), after the optimal adjustment iscomputed from Eq. Eq. (14). The phase differences converge to ¯ x (2)min . (D) The functional pattern associatedwith the phase differences in panel (C). Notice that the only rows and columns that change from the functionalpattern in panel R are the second row and second column. This is due to the fact that only x = θ − θ differs between the two equilibria ¯ x (1)min and ¯ x (2)min . brain regions b r a i n r e g i o n s difference: F − R − m ag n i t ud e Fig 7.
The difference between the functional connectivity F and the functional pattern R . The entries with thelargest magnitude are ≈ .
08, highlighting the stark similarity between the two correlation patterns R and F .23
10 201020 loads l oa d s adjacency matrix A m ag n i t ud e . · B
10 201020 loads l oa d s fault m ag n i t ud e . · C
10 201020 loads l oa d s intervention m ag n i t ud e - · · Fig 8.
Matrices that describe the network interconnections, the fault, and the local intervention of the powernetwork parameters to recover the pre-fault power distribution. (A) The adjacency matrix used in theKuramoto model to simulate the IEEE 39 power network. (B) The fault that disconnects loads 13 and 14. (C)The intervention is localized, in the sense that only branches of the loads connected to the ones affected by thefault and their immediate neighbors require adjustments. The sparsity of the local intervention is promoted bythe usage of the (cid:96)1