Functional control of oscillator networks

Oscillatory activity is ubiquitous in natural and engineered network systems. The interaction scheme underlying interdependent oscillatory components governs the emergence of network-wide patterns of synchrony that regulate and enable complex functions. Yet, understanding, and ultimately harnessing, the structure-function relationship in oscillator networks remains an outstanding challenge of modern science. Here, we address this challenge by presenting a principled method to prescribe exact and robust functional configurations from local network interactions through optimal tuning of the oscillators’ parameters. To quantify the behavioral synchrony between coupled oscillators, we introduce the notion of functional pattern, which encodes the pairwise relationships between the oscillators’ phases. Our procedure is computationally efficient and provably correct, accounts for constrained interaction types, and allows to concurrently assign multiple desired functional patterns. Further, we derive algebraic and graph-theoretic conditions to guarantee the feasibility and stability of target functional patterns. These conditions provide an interpretable mapping between the structural constraints and their functional implications in oscillator networks. As a proof of concept, we apply the proposed method to replicate empirically recorded functional relationships from cortical oscillations in a human brain, and to redistribute the active power flow in different models of electrical grids.


Introduction
Complex coordinated behavior of oscillatory components is linked to the function of many natural and technological network systems [15,17,58].For instance, distinctive network-wide patterns of synchrony [31,71,72] determine the coordinated motion of orbiting particle systems [87], promote successful mating in populations of fireflies [12], regulate the active power flow in electrical grids [23], predict global climate change phenomena [79], dictate the structural development of mother-of-pearl in mollusks [9], and enable numerous cognitive functions in the brain [14,26].Since this rich repertoire of patterns emerges from the properties of the underlying interaction network [73], controlling the collective configuration of interdependent units holds a tremendous potential across science and engineering [75].Despite its practical significance, a comprehensive method to enforce network-wide patterns of synchrony by intervening on the network's structural parameters does not yet exist.
In this work, we develop a rigorous framework that allows us to optimally control the spatial organization of the network components and their oscillation frequencies to achieve desired patterns of synchrony.We abstract the rhythmic activity of a system as the output of a network of diffusively coupled oscillators [4,22] with Kuramoto dynamics.This modeling choice is motivated by the rich dynamical repertoire and wide adoption of Kuramoto oscillators [2].Specifically, we consider an undirected network G = {O, E} of n oscillators with dynamics where ω i ∈ R and θ i ∈ S 1 are the frequency and phase of the i-th oscillator, respectively, A = [A ij ] is the weighted adjacency matrix of G, and O = {1, . . ., n} and E ⊆ O × O denote the oscillator and interconnection sets, respectively.In this work we consider the case where the network G admits both cooperative (i.e., A ij > 0) and competitive (i.e., A ij < 0) [34] interactions among the oscillators, as well as the more constrained case of purely cooperative interactions that arises in several real-word systems.For instance, negative interactions are not physically meaningful in networks of excitatory neurons, in power distribution networks (where the interconnection weight denotes conductance and susceptance of a transmission line), and in urban transportation networks (where interconnections denote the number of vehicles on a road with respect to its maximum capacity).
To quantify the pairwise functional relations between oscillatory units, and inspired by the work in Ref. [5], we define a local correlation metric that, compared to the classical Pearson correlation coefficient, does not 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 such phase differences over time.The right panel illustrates the same oscillator network after a selection of coupling strengths and natural frequencies have been tuned (in red, the structural parameters A and ω are adjusted to A c and ω c ) to obtain the phase differences 2π , which encode the desired functional pattern in the bottom right.In this example, we have computed the closest set (in the 1 -norm sense) of coupling strengths and natural frequencies to the original ones that enforce the emergence of the target pattern.Importantly, only a subset of the original parameters has been modified.depend on sampling time and is more convenient when dealing with periodic phase signals (see Supplementary Information).Given a pair of phase oscillators i and j with phase trajectories θ i (t) and θ j (t), we define the correlation coefficient where < • > t denotes the average over time.A functional pattern is formally defined as the symmetric matrix R whose i, j-th entry equals ρ ij .Importantly, a functional pattern explicitly encodes the pairwise, local, correlations across all of the oscillators, which are more informative than a global observable (e.g., the order parameter [4,6]).
It is easy to see that, if two oscillators i and j synchronize after a certain initial transient, ρ ij → 1 as time increases.
If two oscillators i and j become phase-locked (i.e., their phase difference remains constant over time), then their correlation coefficient converges to some constant value with magnitude smaller than 1.If the phases of two oscillators i and j evolve independently, then their correlation value remains small over time.A few questions arise naturally, which will be answered in this paper.Are all functional patterns achievable?Which network configurations allow for the emergence of multiple target functional patterns?And, if a certain functional pattern can be achieved, is it robust to perturbations? Surprisingly, we reveal that controlling functional patterns can be cast as a convex optimization problem, whose solution can be characterized explicitly.Fig. 1 shows our framework and an example of control of functional pattern for a network with 7 oscillators.In the paper, we will validate our methods by replicating functional patterns from brain recordings in an empirically reconstructed neuronal network, and by controlling the active power distribution in multiple models of power grid.
While synchronization phenomena in oscillator networks have been studied extensively (e.g., see [10,52,57,70,74]), the development of control methods to impose desired synchronous behaviors has only recently attracted the attention of the research community [25,28,45,55].Perhaps the work that is closest to our approach is Ref. [25], where the authors tailor interconnection weights and natural frequencies to achieve a specified level of phase cohesiveness in a network of Kuramoto oscillators.Our work improves considerably upon this latter study, whose goal is limited to prescribing an upper bound to the phase differences, by enabling the prescription of pairwise phase differences and by investigating the stability properties of different functional patterns.Taken together, existing results highlight the importance of controlling distinct configurations of synchrony, but remain mainly focused on the control of "macroscopic" observables of synchrony (e.g., the average synchronization level of all the oscillators).In contrast, our control method prescribes desired pairwise levels of correlation across all of the oscillators, thus enabling a precise "microscopic" description of functional interactions.

Feasible functional patterns in positive networks
A functional pattern is an n×n matrix whose entries are the time-averaged cosine of the differences of the oscillator phases (see equation ( 2)).When the oscillators reach an equilibrium, the differences of the oscillator phases become constant, and the network evolves in a phase-locked configuration.In this case, the functional pattern of the network also becomes constant, and is uniquely characterized by the phase differences at the equilibrium configuration.In this work we study functional patterns for the special case of phase-locked oscillators and, given the equivalence between a desired functional pattern and a desired set of phase differences at equilibrium, convert the problem of generating a functional pattern into the problem of ensuring a desired phase-locked configuration.
We recall that, while convenient for the analysis, phase-locked configurations play a crucial role in the functioning of many natural and man-made networks [38,39,81].
For the undirected network G = (O, E), let x ij = θ j − θ i be the difference of the phases of the oscillators i and j, and let x ∈ R |E| be the vector of all phase differences with (i, j) ∈ E and i < j.The network dynamics (1) can be conveniently rewritten in vector form as (see Methods) where B ∈ R n×|E| is the (oriented) incidence matrix of the network G, D(x) ∈ R |E|×|E| is a diagonal matrix of the sine functions in equation (1), and δ ∈ R |E| is a vector collecting all the weights A ij with i < j.Because we focus on phase-locked trajectories, all oscillators evolve with the same frequency and the vector θ satisfies θ = ω mean 1, where ω mean = 1 n n i=1 ω i is the average of the natural frequencies of the oscillators.Further, since G contains only n oscillators, any phase difference x ij can always be written as a function of n−1 independent differences; for instance, {x 12 , x 23 , . . ., x n−1,n }.In fact, for any pair of oscillators i and j with i < j, it holds x ij = j−1 k=i x k,k+1 .This implies that the vector of all phase differences in equation (3), and in fact any n × n functional pattern, has only n − 1 degrees of freedom and can be uniquely specified with a set of n − 1 independent differences x desired (see Methods).Following this reasoning and to avoid cluttered notation, let ω and notice that the problem of enforcing a desired functional pattern simplifies to (i) converting the desired functional pattern to the corresponding phase differences x desired , and (ii) computing the network weights δ to satisfy the following equation: where we note that the vector ω has zero mean and that, with a slight abuse of notation, D(x) denotes the |E|-dimensional diagonal matrix of the sine of the phase differences uniquely defined by the (n − 1)-dimensional vector x desired .
We begin by studying the problem of attaining a desired functional pattern using only nonnegative weights.
With the above notation, for a desired functional pattern corresponding to the phase differences x, this problem The feasibility of the latter equation, in turn, depends on the projections of the natural frequencies ω on the columns of B: a nonnegative solution exists if ω belongs to the cone generated by the columns of B. This also implies that, if a network admits a desired functional pattern x then, by tuning its weights, the same network can generate any other functional pattern x new such that sign(D(x new )) = sign(D(x)).Thus, by properly tuning its weights, a network can generally generate a continuum of functional patterns determined uniquely by signs of its incidence matrix and the oscillators natural frequencies.This property is illustrated in Fig. 2 for the case of a line network.
A sufficient condition for the feasibility of Problem ( 5) is as follows: There exists δ ≥ 0 such that BD(x)δ = ω if there exists a set S satisfying: (i.a) D ii (x)D jj (x)B T :,i B :,j ≤ 0 for all i, j ∈ S with i = j and D ii , D jj = 0; (i.b) ω T B :,i D ii (x) > 0 for all i ∈ S; Equivalently, the above conditions ensure that ω is contained within the cone generated by the columns of The existence of such an Hamiltonian path, together with a positive projection of ω onto B:,H , also ensure that there exists a strictly positive δ > 0 solution to BD(x)δ = ω.In particular, for any choice of x 12 , x 23 ∈ (0, π), equation ( 4) reveals that if 0 < A 13 < 0.5/ sin(x 12 + x 23 ), then there exist strictly positive weights A 12 > 0 and A 23 > 0 such that δ > 0.
B:,S (see Fig. 3a for a self-contained example).To see this, rewrite the pattern assignment problem BD(x)δ = ω as where the subscripts S and S denote the entries corresponding to the set S and the remaining ones, respectively.
If the vectors B :,i , i ∈ S, are linearly independent, condition (i.a) implies that D S,S B T :,S B :,S D S,S is an M-matrix; that is, a matrix which has nonpositive off-diagonal elements and positive principal minors [36].Otherwise, the argument holds verbatim by replacing S with any subset S ind ⊂ S such that the vectors B :,i , i ∈ S ind , are linearly independent.Condition (i.c) guarantees the existence of a solution to B :,S D S,S δ S = ω.A particular solution to the latter equation is δ S = (B :,S D S,S (x)) † ω = (D S,S (x)B T :,S B :,S D S,S (x)) −1 D S,S (x)B T S,S ω > 0 where (•) † denotes the Moore-Penrose pseudo-inverse of a matrix.The positivity of δ S follows from condition (i.b) and the fact that the inverse of an M-matrix is element-wise nonnegative [36].We conclude that a solution to equation ( 6) is given by δ S = (B :,S D S,S (x)) † ω > 0 and δ S = 0.
To avoid disconnecting edges or to maintain a fixed network topology, a functional pattern should be realized in Problem (5) using a strictly positive weight vector (that is, δ > 0 rather than δ ≥ 0).While, in general, this is a considerably harder problem, a sufficient condition for the existence of a strictly positive solution δ > 0 is that the network with incidence matrix B contains an Hamiltonian path, that is, a directed path that visits all the oscillators exactly once (Fig. 3b shows a network containing an Hamiltonian path).Namely, There exists a strictly positive solution δ > 0 to BD(x)δ = ω if (ii.a) the network with incidence matrix B contains a directed Hamiltonian path H; The incidence matrix B:,H of a directed Hamiltonian path H has two key properties.First, it comprises n − 1 linearly independent columns, since the path covers all the vertices and contains no cycles.This guarantees that ω ∈ Im(B :,H ).Second, the columns of the incidence matrix B:,H satisfy BT :,i B:,i = 2 and BT :,i B:,j ∈ {0, −1} for all i, j ∈ H, i = j.Then, letting the set S in the result above identify the columns of the Hamiltonian path, conditions (ii.a) and (ii.b) imply (i.a)-(i.c),thus ensuring the existence of a nonnegative set of weights δ that solves the pattern assignment problem BD(x)δ = ω.Furthermore, by rewriting the pattern assignment problem as in equation ( 6), the following vector of interconnection weights solves such equation (see Methods): † ω − B :, HD H, H(x)δ H .
Because B:,H = B :,H D H,H defines an Hamiltonian path and because of (ii.b), the vector (B :,H D H,H (x)) † ω contains only strictly positive entries.Thus, for any sufficiently small and positive vector δ H, the weights δ H are also strictly positive, ultimately proving the existence of a strictly positive solution to the pattern assignment problem (see Methods).Fig. 3c illustrates a self-contained example.
Taken together, the results presented in this section reveal that the interplay between the network structure and the oscillators' natural frequencies dictates whether a desired functional pattern is achievable under the constraint of nonnegative (or even strictly positive) interconnections.First, dense positive networks with a large number of edges are more likely to generate a desired functional pattern, since their incidence matrix features a larger number of candidate vectors to satisfy conditions (i.a)-(i.c).Second, densely connected networks are also more likely to contain an Hamiltonian path, thus promoting also strictly positive network designs.Third, after an appropriate relabeling of the oscillators such that any interconnection from i to j in the Hamiltonian path satisfies i < j, condition (ii.b) is equivalent to requiring that ω i < ω j .That is, the feasibility of a functional pattern is guaranteed when the natural frequencies increase monotonically with the ordering identified by the Hamiltonian path.This also implies, for instance, that sparsely connected positive networks, and not only dense ones, can attain a large variety of functional patterns.An example is a connected line network with increasing natural frequencies, which can generate, among others, any functional pattern defined by phase differences that are smaller than π 2 (trivially, when the phase differences are smaller than π 2 and the natural frequencies are increasing, a line network contains an Hamiltonian path and the vector of natural frequencies has positive projections onto the columns of the incidence matrix).Fig. 2a contains an example of such a network.

Compatibility of multiple functional patterns
A single choice of the interconnections weights can allow for multiple desired functional patterns, as long as they are compatible with the network dynamics in equation (1).In this section, we provide a characterization of compatible functional patterns in a given network, and derive conditions for the existence of a set of interconnection weights that achieve multiple desired functional patterns.Being able to concurrently assign multiple functional patterns is crucial, for instance, to the investigation and design of memory systems [68], where different patterns of activity correspond to distinct memories.Furthermore, our results complement previous work on the search of equilibria in oscillator networks [44].
To find a set of functional patterns that exists concurrently in a given network with fixed interconnection weights δ, we exploit the algebraic core of equation ( 4) and show that the kernel of the incidence matrix B uniquely determines the equilibria of the network.In fact, for a given network (i.e., δ with nonzero components) all compatible equilibria x (i) , i ∈ {1, . . ., }, must satisfy From equation (7), we can see that the sine vector of all compatible equilibria must belong to a specific affine subspace of R |E| : Rewriting equation (4) in the above form connects the existence of distinct functional patterns with ker(B), the latter featuring a number of well known properties.For instance, it holds that dim(ker(B)) = |E| − n + c, where c The fourth panel reveals that the two functional patterns compatible with x (j) , j ∈ {1, . . ., 4}, correspond to x (5) = [2π/3 2π/3] T and x (6) = [−2π/3 − 2π/3] T (in red).
is the number of connected components in a network, and that ker(B) coincides with the subspace spanned by the signed path vectors of all undirected cycles in the network [30].Notice also that, after a suitable reordering of the phase differences in x, we can write sin(x) = [sin(x T desired ) sin(x T dep )] T , where x dep denotes the phase differences dependent on n − 1 desired phase differences x desired .Thus, all the x desired for which sin(x dep ) intersects the affine space described by diag(δ) −1 B † ω + ker(B) identify compatible functional patterns.
To showcase how the intimate relationship between the network structure and the kernel of its incidence matrix enables the characterization of which (and how many) compatible patterns coexist, we consider three essential network topologies: trees, cycles, and complete graphs.For the sake of simplicity, we let δ = 1 and ω = 0, so that equation (8) holds whenever sin(x (i) ) ∈ ker(B).In networks with tree topologies it holds ker(B) = 0, and sin(x (i) ) = 0 is satisfied by 2 n−1 patterns of the form x (i) jk = 0, π, for all (j, k) ∈ E. Consider now cycle networks, where ker(B) = span 1.For any cycle of n ≥ 3 oscillators, two families of patterns are straightforward to derive.First, there are 2 n−1 patterns of the form x (i) k,k+1 = 0, π, with k = 1, . . ., n − 1, and k,k+1 .Second, there are n − 1 splay states [22], where the oscillators' phases evenly span the unitary circle, with x n , m = 1, . . ., n − 1, (j, k) ∈ E. Fig. 4 illustrates the compatible functional patterns satisfying equation (8) in a positive network of three fully synchronizing oscillators.In general, however, cycle networks of identical oscillators admit infinite coexisting patterns.For instance, Fig. 5 shows how we can parameterize infinite equilibria with a scalar γ ∈ S 1 in a cycle of n = 4 oscillators.Finally, as complete graphs . . .A homogeneous cycle network admits infinite compatible functional patterns.Since ker(B) = span 1, the cycle network admits infinite compatible equilibria, which can be parameterized by γ ∈ S 1 as Any arbitrarily small variation of γ yields sin(x (i) (γ)) ∈ ker(B).The right panel illustrates the patterns associated with x (i) (γ), i = 1, . . ., 5 for increments of γ of 0.2 radians.are equivalent to a composition of cycles, they also admit infinite compatible patterns that can be parameterized akin to what occurs in a simple cycle (see Supplementary Figure 11).
We now turn our attention to finding the interconnection weights that simultaneously enable a collection of ≥ 1 desired functional patterns {x (i) } i=1 .We first notice that equation (7) reveals that to achieve a desired functional pattern x (i) with components not equal to kπ, k ∈ Z, the network weights δ must belong to an |E|-dimensional affine subspace of R |E| : Let Γ i = D(x (i) ) −1 B † ω + ker(B) .Then, to concurrently realize a collection of patterns {x (i) } i=1 , a solution to equations ( 9) exists if and only if i=1 Γ i = ∅.It is worth noting that, whenever the latter intersection coincides with a singleton, then there exists a single choice of network weights that realizes {x (i) } i=1 .However, if i=1 Γ i corresponds to a subspace, then infinite networks can realize the desired collection of functional patterns.We conclude by emphasizing that a positive δ that achieves the desired patterns exists if and only That is, if the network weights belong to the nonempty intersection of the affine subspaces with the positive orthant.

Stability of functional patterns
A functional pattern is stable when small deviations of the oscillators phases from the desired configuration lead to vanishing functional perturbations.Stability is a desired property since it guarantees that the desired functional pattern is robust against perturbations to the oscillators dynamics.To study the stability of a functional pattern, we analyze the Jacobian of the Kuramoto dynamics at the desired functional configuration, which reads as [22] where L(x) denotes the Laplacian matrix of the network with weights scaled by the cosines of the phase differences (the weight between nodes i and j is A ij cos(θ j − θ i )).The functional pattern x is stable when the eigenvalues of the above Jacobian matrix have negative real parts.For instance, if all phase differences are strictly smaller than π 2 (that is, the infinity-norm of x satisfies x ∞ < π 2 ), then the Jacobian in equation ( 10) is known to be stable [22].In the case that both cooperative and competitive interactions are allowed, we can ensure stability of a desired pattern by specifying the network weights in δ such that 2 and A ij < 0 otherwise, so that the matrix L becomes the Laplacian of a positive network (see Methods).Furthermore, we observe that in the particular case where some differences |x ij | = π 2 , the network may become disconnected since cos(x ij ) = 0.
Because the Laplacian of a disconnected network has multiple eigenvalues at zero, marginal stability may occur, and phase trajectories may not converge to the desired pattern.
When some phase differences are larger than π 2 and the network allows only for nonnegative weights, then stability of a functional pattern is more difficult to assess because the Jacobian matrix becomes a signed Laplacian [86].The off-diagonal entries of a signed Laplacian satisfy L ij > 0 whenever |x ij | > π 2 , thus possibly changing the sign of its diagonal entries L ii = j A ij cos(x ij ) and violating the conditions for the use of classic results from algebraic graph theory for the stability of Laplacian matrices.To derive a condition for the instability of the Jacobian in equation ( 10), we exploit the notion of structural balance.We say that the cosine-scaled network with Laplacian matrix L is structurally balanced if and only if its oscillators can be partitioned into two sets, O 1 and O 2 , such that all (i, j) ∈ E with A ij cos(x ij ) < 0 connect oscillators in O 1 to oscillators in O 2 , and all then its Laplacian has mixed eigenvalues [86].Therefore, we conclude the following: If the functional pattern x yields a structurally balanced cosine-scaled network, then x is unstable.
The above condition allows us to immediately assess the instability of functional patterns for the special cases of line and cycle networks.In fact, for a line network with positive weights, x is unstable whenever |x ij | > π 2 for any i, j.Instead, for a cycle network with positive weights, the pattern x can be stable only if it contains at most one phase difference π 2 < |x ij | < γ, where γ ≈ 1.789776 solves γ − tan(γ) = 2π (see Supplementary Information).In the next section, we propose a heuristic procedure to correct the interconnection weights in positive networks to promote stability of a functional pattern.

Optimal interventions for desired functional patterns
Armed with conditions to guarantee the existence of positive interconnections that enable a desired functional pattern, we now show that the problem of adjusting the network weights to generate a desired functional pattern can be cast as a convex optimization problem.Formally, for a desired functional pattern x and network weights δ, we seek to solve min α α 2 (11) subject to where α ∈ R |E| are the controllable modifications of the network weights, and • 2 denotes the 2 -norm.Fig. 6a illustrates the control of a functional pattern in a line network of n = 4 oscillators.
The minimization problem ( 11) is convex, thus efficiently solvable even for large networks, and may admit multiple minimizers, thus showing that different networks may exhibit the same functional pattern.Moreover, in light of our results above, Problem ( 11) can be easily adapted to assign a collection of desired patterns {x (i) } i=1 .
To do so, we simply replace the constraint (11a) with Fig. 6b illustrates an example where we jointly allocate two functional patterns for a complete graph with n = 7 oscillators (see Supplementary Information for more details on this example).
Note that the minimization problem (11) does not guarantee that the functional pattern x is stable for the network with weights δ + α * .To promote stability of the pattern x, we use a heuristic procedure based on the classic Gerschgorin's theorem [29].Recall that stability of x is guaranteed when the associated Jacobian matrix has a Laplacian structure, with negative diagonal entries and nonnegative off-diagonal entries.Further, instability of x depends primarily on the negative off-diagonal entries A ij cos(x ij ) of the Jacobian (these entries T is associated to interconnection weights δ = [3.40263.4641 6.4721] T .Applying the optimal correction α * yields positive interconnection weights δ + α * = [6.47213.4026 3.4641] T that achieve the desired functional patterns x desired .b Joint allocation of two compatible equilibria for the phase difference dynamics.By taking θ 1 as a reference, we choose two points for the phase differences x 1i = θ i − θ 1 , i ∈ {2, . . ., 7}, to be imposed as equilibria in a network of n = 7 oscillators: x T .In this example, we find a set of interconnection weights (δ + α * ) that solves the minimization problem (11) with constraint (12).The trajectories start at the (unstable) equilibrium point x (1) desired at time t = 0, and converge to the point x (2) desired after a small perturbation p ∈ T 7 , with p i ∈ [0 0.05], is applied to the phase difference trajectories at time t = 50.are negative because the sign of the network weight A ij is different from the sign of the cosine of the desired phase difference x ij ).Thus, reducing the magnitude of such entries A ij heuristically moves the eigenvalues of the Jacobian towards the stable half of the complex plane (this phenomenon can be captured using the Gerschgorin circles, as we show in Fig. 7 for a network with 7 nodes).To formalize this procedure, let δ N and α N denote the entries of the weights δ and tuning vector α, respectively, that are associated to negative interconnections A ij cos(x ij ) < 0 in the cosine-scaled network.Then, the optimization problem that enacts the proposed strategy becomes: and (δ + α) ≥ 0, Carefully reducing the weights δ N + α N promotes stability of the target pattern.Fig. 7 illustrates the shift of the Jacobian's eigenvalues while the optimal tuning vector α * is gradually applied to a 7-oscillator network to achieve stability of a functional pattern containing negative correlations (the network parameters can be found T , where The left plot illustrates the Gerschgorin disks (in blue) and the Jacobian's eigenvalues locations for the original network (as dark dots).The complex axis is highlighted in purple.It can be observed in the zoomed-in panel that one eigenvalue is unstable (λ 2 = 0.0565, in red).The optimal correction α * is gradually applied to the existing interconnections from the left-most panel to the right-most one at 1   3   increments.The right zoomed-in panel shows that, as a result of our procedure, n − 1 eigenvalues ultimately lie in the left-hand side of the complex plane (λ 1 = 0 due to rotational symmetry and λ 2 = −0.0178, in green). in the Supplementary Information).Finally, we remark that the procedure in ( 21) can be further refined by introducing scaling constants to penalize δ N + α N 2 differently from the modification of other interconnection weights (see Supplementary Information for further details and an example).
The minimization problems (11) and (21) do not allow us to tune the oscillators' natural frequencies, and are constrained to networks with positive weights.When any parameter of the network is unconstrained and can be adjusted to enforce a desired functional pattern, the network optimization problem can be generalized as min where β denotes the correction to the natural frequencies.Problem ( 14) always admits a solution because β can be chosen to satisfy the constraint (25a) for any choice of the network parameters δ + α.Further, the (unique) solution to the minimization problem (14) can also be computed in closed form: where I n denotes the n × n identity matrix.
We conclude this section by noting that the minimization problems ( 11)-( 14) can be readily extended to include other vector norms besides the 2 -norm in the cost function (e.g., the 1 -norm to promote sparsity of the corrections), and to be applicable to directed networks.The latter extension can be attained by replacing the constraints (11a) and (25a) with a suitable rewriting of the matrix form (4). We refer the interested reader to the Supplementary Information for a comprehensive treatment of this extension and an example.

Applications to complex networks
In the remainder of this paper we apply our methods to an empirically reconstructed brain network and to the IEEE 39 New England power distribution network.In the former case, we use the Kuramoto model to map structure to function, and find that local metabolic changes underlie the emergence of functional patterns of recorded neural activity.In the latter case, we use our methods to restore the nominal network power flow after a fault.

Local metabolic changes govern the emergence of distinct functional patterns in the brain
The brain can be studied as a network system in which Kuramoto oscillators approximate the rhythmic activity of different brain regions [14,45,59,60].Over short time frames, the brain is capable of exhibiting a rich repertoire 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 as biomarkers for different psychiatric and neurological disorders [65].
To shed light on the structure-function relationship of the human brain, we utilize Kuramoto oscillators evolving on an empirically reconstructed brain network.We hypothesize that the intermittent emergence of diverse patterns stems from changes in the oscillators' natural frequencies -which can be thought of as endogenous changes in metabolic regional activity regulated by glial cells [27] or exogenous interventions to modify undesired synchronization patterns [45].First, we show that phase-locked trajectories of the Kuramoto model in equation (1) can be accurately extracted from noisy measurements of neural activity and are a relatively accurate approximation of neural activity.
We employ structural (i.e., interconnections between brain regions) and functional (i.e., time series of recorded neural activity) data from Ref. [60].Structural connectivity consists of a sparse weighted matrix A whose entries represent the strength of the physical interconnection between two brain regions.Functional data comprise time series of neural activity recorded through functional magnetic resonance imaging (fMRI) of healthy subjects as they rest.Because the phases of the measured activity have been shown to carry most of the information contained in the slow oscillations recorded through fMRI time series, we follow the steps in Ref. [60] to obtain such phases from the data by filtering the time series in the [0.04, 0.07] Hz frequency range (Supplementary Information).Next, since frequency synchronization is thought to be a crucial enabler of information exchange between different brain regions and homeostasis of brain dynamics [53,82], we focus on functional patterns that arise from phase-locked trajectories, as compatible with our analysis.For simplicity, we restrict our study to the cingulo-opercular cognitive system, which, at the spatial scale of our data, comprises n = 12 interacting brain regions [61].
We identify time windows in the filtered fMRI time series where the signals are phase-locked, and compute two matrices for each time window: a matrix F ∈ R 12×12 of Pearson correlation coefficients (also known as functional connectivity), and a functional pattern R ∈ R 12×12 (as in equation ( 2)) from the phases extracted by solving the nonconvex phase synchronization problem [11].Strikingly, we find that F −R 2 1 consistently (see Supplementary Information and Supplementary Figure 16b), thus demonstrating that our definition of functional pattern is a good replacement of the classical Pearson correlation arrangements in networks with oscillating states.Building upon the above finding, we test whether the oscillators' natural frequencies embody the parameter that links the brain network structure to its function (i.e., structural and functional matrices).We set ω = We obtain R − F 2 = 0.2879.Additionally, we verify that the analysis of the Jacobian spectrum (see equation ( 10)) accurately predicts the stability of the phase-locked trajectories.Supplementary Figure 7a illustrates the basin of attraction of R, which we numerically estimate to be half of the torus.
BD(x)δ, where x are phase differences obtained from the previous step, and integrate the Kuramoto model in equation (1) with random initial conditions close to x.We show in Fig. 8 that the assignment of natural frequencies according to the extracted phase differences promotes spontaneous synchronization to accurately replicate the empirical functional connectivity F .
These results corroborate the postulate that structural connections in the brain support the intermittent activation of specific functional patterns during rest through regional metabolic changes.Furthermore, we show that the Kuramoto model represents an accurate and interpretable mapping between the brain anatomical organization and the functional patterns of frequency-synchronized neural co-fluctuations.Once a good mapping is inferred, it can be used to define a generative brain model to replicate in silico how the brain efficiently enacts large-scale integration of information, and to develop personalized intervention schemes for neurological disorders related to abnormal synchronization phenomena [47,83].

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 infrastructure.Modern, reconfigurable power networks are expected to be resilient to distributed faults and malicious cyber-physical attacks [84] while being able to rapidly adapt to varying load demands.In addition, climate change is straining service reliability, as underscored by recent events such as the Texas grid collapse after Winter Storm Uri in February 2021 [13], and the New Orleans blackout after Hurricane Ida in August 2021 [1].Therefore, there exists a dire need to design control methods to efficiently operate these networks and react to unforeseen disruptive events.
The Kuramoto model in equation (1) has been shown to be particularly relevant in the modeling of large The fault causes the voltage phases θ to depart from normal operating conditions, which could 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 and promote a local (sparse) intervention, we solve the optimization in equation ( 11) by minimizing the 1-norm of the structural parameter modification δ with no scaling parameters in the cost functional.The network returns to the initial operative conditions with a localized modification of the neighboring transmission lines' impedances.
distribution networks and microgrids [23].Preliminary work on the control of frequency synchronization in electrical grids modeled through Kuramoto oscillators has been developed in Ref. [67].Here, we present a method that leverages our findings to guarantee not only frequency synchronization, but also a target pattern of active power flow.Our method can be used for power (re)distribution with respect to specific pricing strategies, fault prevention (e.g., when a line overheats) and recovery (e.g., after the disconnection of a branch).Furthermore, thanks to the formal guarantees that our method prescribes, we are able to prevent Braess' Paradox in power networks [85], which is a phenomenon where the addition of interconnections to a network may impede its synchronization.
It has been shown in Ref. [23,Lemma 1] that the load dynamics (nodes 1-29 in Fig. 9a) in a structurepreserving power grid model have the same stable synchronization manifold of equation (1).In this model, where D i = 1, which is possibly due to local excitation controllers.Notice that, when the phase angles θ are phase-locked and A ij is fixed, the active power flow is given by A ij sin(θ j − θ i ).
We posit that solving the problem in equation (11) to design a local reconfiguration of the network parameters can recover the power distribution before a line fault or provide the smallest parameter changes to steer the load powers to desired values.In practice, control devices such as Flexible Alternating Current Transmission Systems (FACTs) allow operators and engineers to change the network parameters [42,50].We demonstrate the effectiveness of our approach by recovering a desired power distribution in the IEEE 39 New England power distribution network after a fault.During a regime of normal operation, we simulate a fault by disconnecting the line between two loads and solve the problem in equation ( 11) to find the minimum modification of the couplings a ij that recovers the nominal power distribution.
We first utilize MATPOWER [89] to solve the power flow problem.Then we use the active powers p 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 equation (1).
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 pattern across the loads.Next, to simulate a line fault, we disconnect one line.By solving the problem in equation ( 11) with x desired corresponding to the pre-fault steady-state voltage phase differences, we compute the smallest variation of the remaining parameters (i.e., admittances) so that the original functional pattern can be recovered.Fig. 9b illustrates the effectiveness of our procedure at recovering the nominal pattern of active power flow by means of minimal and localized interventions (see also Supplementary Figure 17).
The above application is based on a classical lossless structure-preserving power network model [23].However, in the power systems literature, more complex dynamics that relax some of the modeling assumptions have been proposed.For instance, Ref. [66] uses a third-order model (or one-axis model) that takes into account transient voltage dynamics.Instead, Ref. [33] studies the case of interconnections with power losses, which lead to a network of phase-lagged Kuramoto-Sakaguchi oscillators [63].We show in the Supplementary Information that our procedure to recover a target functional pattern can still be applied successfully to a wide range of situations involving these more realistic models.

Contribution and future directions
Distinct configurations of synchrony govern the functioning of oscillatory network systems.This work presents a simple and mathematically grounded mapping between the structural parameters of arbitrary oscillator networks and their components' functional interactions.
The tantalizing idea of prescribing patterns in networks of oscillators has been investigated before, yet only partial results had been reported in the literature.Here, we demonstrate that the control of patterns of synchrony can be cast as optimal (convex) design and tuning problems.We also investigate the feasibility of such optimizations in the cases of networks that admit negative coupling weights and networks that are constrained to positive couplings.Our control framework also allows us to prescribe multiple desired equilibria in Kuramoto networks, a problem that is relevant in practice and had not been investigated before.
As stability of a functional pattern may be a compelling property in many applications, we explore conditions to test and enforce the stability of functional patterns.We show that such conditions are rather straightforward in the case of networks that admit both cooperative and competitive interactions.However, stability of target functional patterns in networks that are constrained to only cooperative interactions is a more challenging task, for which we demonstrate that any pattern associated to a structurally-balanced cosine-weighted network cannot be stable.To overcome this issue, we propose an heuristic procedure that adjusts the oscillators' coupling strengths to violate the structurally-balanced property and promote the stability of functional patterns with negative correlations.While heuristic, our procedure for stability has proven successful in all our numerical studies.
Notice that, differently from methods that study an "average" description of the system at near-synchronous state (see, e.g., Ref. [76]), here we assess the stability of exact target phase-locked trajectories where phase values can instead be arbitrarily spread over the torus.This method can also be extended to assess the stability of equilibria of higher-dimensional oscillators, provided that the considered equilibrium state is a fixed point.
We emphasize that our results are also intimately related to the long standing economic problem of enhancing network operations while optimizing wiring costs.In any complex system where synchrony between components ensures appropriate functions, it is beneficial to maximize synchronization while minimizing the physical variations of the interconnection weights [51].Compatible with this principle, neural systems are thought to have evolved to maximize information processing by promoting synchronization through optimal spatial organization [3].Inspired by the efficiency observed in neurobiological circuitry, equation ( 11) could be utilized for the design optimal interaction schemes in large-scale computer networks whose performance relies on synchronizationbased tasks [41].
An important consideration that highlights the general contributions of the present study is that being able to specify pairwise functional relationships between the oscillators also solves problems such as phase-locking, full, and cluster synchronization.Yet, the converse is not true.In fact, even in the general setting of cluster synchronization [56] -where distinct groups of oscillators behave cohesively -one can only achieve a desired synchronization level within the same cluster, but not across clusters, which instead is possible with our approach.
Specifically, in cluster-synchronization oscillators belonging to the same cluster are forced to synchronize, thus implying that the associated diagonal blocks of the functional pattern R display values close to 1. Seminal work in Ref. [40] developed a nonlinear feedback control to change the coupling functions in equation ( 2) to engineer clusters of synchronized oscillators, whereas the authors of Ref. [24] propose the formation of clusters through selective addition of interconnections to the network.More recently, the control of partially synchronized states with applications to brain networks is studied in Ref. [45] by means of structural interventions, and in Ref. [8] via exogenous stimulation.Ultimately, the results presented in this work not only complement but go well beyond the control of macroscopic synchronization observables and partial synchronization by allowing to specify the synchrony level of pairwise interactions.
While our contributions include the analysis of the stability properties of functional patterns, here we do not assess their basins of attraction.We emphasize that, in general, the estimation of the basin of attraction of the attractors of nonlinear systems remains an outstanding problem, and even the most recent results rely on numerical approaches or heavy modeling assumptions [16].Further, in the case of coupled Kuramoto oscillators, existing literature shows that the number of equilibria for the phase differences increases significantly with the cardinality of the network [44], making the study of basins of attraction extremely challenging.In the Supplementary Information, we extend previous work on identical oscillators to networks with heterogeneous oscillators, and show that that functional patterns can be at most almost-globally stable in cluster-synchronized positive networks.Yet, a precise estimation of the basin of attraction for any arbitrary target pattern goes beyond the scope of this work and may require the derivation of ad-hoc principles based on Lyapunov's stability theory [88].
The framework presented in this work has other limitations, which can be addressed in follow up studies.
First, despite its capabilities in modeling numerous oscillatory network systems, the Kuramoto model cannot capture the amplitude of the oscillations, making it most suitable for oscillator systems where most of the information is conveyed by phase interaction as demonstrated in Ref. [60] for resting brain activity.To model brain activity during cognitively demanding tasks such as learning, higher-dimensional oscillators may be more suitable [62].Second, the use of phase-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 which only a subset of the oscillators is frequency-synchronized.Third and finally, in some situations the network parameters are not fully known.While being still an active area of research, network identification of oscillator systems may be employed in such scenarios [54].
Directions of future research can be both of a theoretical and practical nature.For instance, follow up studies can focus on the derivation of a general condition for the stability of a feasible functional pattern in positive networks.Further, a thorough investigation of which network structures allow for multiple prescribed equilibria may be particularly relevant in the context of memory systems, where different patterns are associated to different memory states.Specific practical applications may also require the inclusion of sparsity constraints on the accessible structural parameters for the implementation of the proposed control and design framework.

Matrix form of equation (1) 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).In particular, i points to j if i < j, and is the source node of the interconnection , B k = 1 if oscillator k is the sink node of the interconnection , and B k = 0 otherwise.The matrix form of equation ( 1) can be written as where D(x) is the diagonal matrix of the sine functions in equation (1).
When the oscillators evolve in a phase-locked configuration, the oscillator frequencies become equal to each other and constant.In particular, since 1 T B = 0, we have :,H ω is strictly positive.Then, the vector (DH,H(x)B T :,H B:,HDH,H(x)) −1 DH,H(x)B T :,H ω is also strictly positive, and so is the solution vector δH for any sufficiently small and positive vector δ H.

Enforcing stability of functional patterns in networks with cooperative and competitive interactions
To ensure that the Jacobian matrix in equation ( 10) is the Laplacian of a positive network and, thus, stable, we solve the problem in equation ( 5) with a slight modification of the constraints.Specifically, we post-multiply the matrix B in equation ( 5) as B∆, where ∆ = diag({sign(cos(xij))} (i,j)∈E ) is a matrix that changes the sign of the columns of B associated to negative weights in the cosine-scaled network: Solving for positive interconnection weights the problem in equation ( 5) under the above modified constraint yields a stable Jacobian in a network where the final couplings are ∆δ.

Heuristic procedure to promote stability of functional patterns in positive networks
We provide a heuristic procedure to promote the stability of functional patterns that include negative correlations in a network with nonnegative weights.Our procedure relies on the definition of Gerschgorin disks and the Gerschgorin Theorem.
Definition of Gerschgorin disk.Let M ∈ C n×n be a complex matrix.The i-th Gerschgorin disk is Di = (Mii, ri), i = 1, . . ., n, where the radius is ri = j =i |Mij| and the center is Mii.
Theorem 2 (Gerschgorin [29]).The eigenvalues of the matrix M lie within the union n i=1 Di of its Gerschgorin disks.
Whenever all target phase differences in x satisfy |xij| ≤ π 2 , the Gerschgorin disks of the Jacobian in equation (10) all lie in the closed left half-plane.However, for patterns x containing phase differences |xij| ≥ π 2 , the union of the Gerschgorin disks intersects the right half-plane.Reducing the magnitude of the entries satisfying Aij cos(xij) < 0 effectively shrinks the radius of the Gerschgorin disks that overlap with the right half-plane and shifts their centers towards the left-half plane due to the structure of the Jacobian matrix.We remark that the procedure in equation ( 21) is a heuristic, and it is provably effective only when all interconnections with Aij cos(xij) < 0 can be removed, so that all the Gerschgorin disks lie completely in the left-half plane.
is stable if and only if for i = 2, . . ., n, with n − i + 3 1 if i = 2, and n + 1 1; and n → ∞, the largest possible value for γ such that x cycle is stable tends to the value γ ≈ 1.789776, which is the solution to γ − tan(γ) = 2π.

Proof of Theorem 2:
To assess the stability of the equilibrium xcycle = [x 12 x 23 . . .
, 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. [86,Corollary IV.7], we have that a necessary and sufficient condition for the Laplacian matrix L(x cycle ) of the cosine-scaled network to be positive semidefinite is with R 12 being the effective resistance of the graph in which the edge (1, 2) has been removed.That is, Since the adjacency matrix satisfies A = A T and x cycle is an equilibrium for the difference dynamics of the cycle network, the network weights must be identical pairwise: for i = 2, . . ., n with the convention n − i + 3 1 if i = 2. Thus, the angles being identical pairwise implies that the network weights must also be identical pairwise, which yields condition (i) of Theorem 2.Moreover, plugging the network weights from equation ( 17) into equation ( 16) yields which makes the condition in Eq. ( 15) become, after algebraic calculations, condition (ii) of Theorem 2: Thus, given that the Laplacian L(x cycle ) is positive definite, the Jacobian J(x cycle ) is stable, and its only zero eigenvalue is due to rotational symmetry of the right-hand side of the Kuramoto dynamics (equation ( 2) of the main manuscript). For Hence, the right-hand side of equation ( 19) becomes cot(ϕ) n−1 , and lim n→∞ cot(ϕ) n−1 = 1 2π−γ .Since |γ| > π 2 , plugging the limit value for 1 2π−γ into equation (19) and solving for the equality yields γ − tan(γ) = 2π, whose unique solution is γ ≈ 1.789776.This concludes the proof.acksquare

Functional patterns possess at most almost-global stability in positive cluster-synchronized networks
In general, the analysis and estimation of the basin of attraction of nonlinear systems remains an outstanding problem, and even the most recent results rely on numerical approaches or heavy modeling assumptions [16].To the best of our knowledge, Ref. [88] provides the most up-to-date study of the basins of attraction of synchronized Kuramoto oscillators.However, the authors in Ref. [88] only derive estimates of the basin of attraction for the fully synchronized case in networks of identical oscillators.Below, we show that functional patterns associated with phase-locked trajectories in cluster-synchronized positive networks are at most almost-globally stable.Our results extend previous work on identical oscillators to the case of oscillators with heterogeneous natural frequencies.
Existing literature shows that the number of equilibria for the phase differences of heterogeneous oscillators evolving on random graphs increases significantly with the cardinality of the network [44].Therefore, we begin by restricting our analysis to the case of identical oscillators.In general, for any connected network of identical oscillators with cooperative connections A ij ≥ 0, there is no unique pattern.In fact, the largest basin of attraction is achieved by S1 -synchronizing graphs (e.g., complete graphs, acyclic graphs, and sufficiently dense graphs [22,78]), which feature almost-global stability of the fully synchronized functional pattern.In this class of networks, the only stable equilibrium x satisfies x ij = 0 for all i, j, and there exists a finite number of unstable equilibria.For instance, consider two connected identical oscillators: ω 1 = ω 2 .The only stable equilibrium for the phase differences dynamics is x 12 = θ 2 − θ 1 = 0. Yet, there also exists an unstable equilibrium x 12 = π.
Note that, in general, almost-global stability cannot be guaranteed for arbitrary classes of sparse networks of identical oscillators, as other synchronization manifolds besides the fully synchronized one can be stable -for example, splay states emerge in Cayley graphs [43].This latter class of graphs admits two stable equilibria (full synchronization and splay states), which implies the coexistence of two stable patterns.
We can generalize the above observations to networks of heterogeneous oscillators.To do so, we leverage cluster synchronization -a phenomenon where distinct groups of synchronized oscillators coexist in a network [46] -of S 1 -synchronizing clusters with possibly different natural frequencies.We consider cluster synchronization  15).2Thus, for a given number of clusters m ≥ 1, each of the n m possible ways to partition the network yields a functional pattern with diagonal blocks corresponding to synchronized clusters.By extending the above observations on the stability of S 1 -synchronizing graphs, we find that, in networks where at least one cluster satisfies |C k | ≥ 2, at most almost-global stability of cluster-synchronized functional patterns can be achieved.In fact, for any choice of intra-cluster phase differences from the finite set of unstable equilibria in S 1 -synchronizing graphs (i.e., x ij = π for at least one intra-cluster phase difference), there exists a value for the inter-cluster phase differences for which the network admits phase-locked trajectories where intra-cluster phase differences satisfy x ij = 0 or x ij = π.

The network parameters read
Note that the two clusters are S 1 -synchronizing graphs.It can be shown that there exists an unstable equilibrium at x desired = [x 12 x 23 x 34 ] T = [π 0.25268 π] T .Thus, whenever the pattern R associated to the partition C is stable, it is at most almost-globally stable.

Allocation of multiple phase-locked equilibria by tailoring of the network structural parameters
The convex optimization problem proposed in the main text allows to tailor the network weights and the natural frequencies of the oscillators to specify multiple equilibria for the dynamics of the phase differences x.These equilibria 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 equilibria for the phase difference dynamics.Specifically, by taking θ 1 as a reference, we choose two points for the phase differences x 1i = θ i − θ 1 to be set as equilibria: x T .
The initial network parameters (adjacency matrix and zero-mean natural frequencies) read as: , respectively.To impose the desired equilibria, we numerically solve the following convex problem through standard cvx routines [32]: The solution α * yields the following corrected adjacency matrix: An investigation of the Jacobian spectrum (see main text for results on stability) of the phase differences computed at the two equilibria x desired and x desired reveals that the first equilibrium point is unstable and that the second one is locally stable.We illustrate the outcome of the above procedure to specify multiple equilibria in Fig. 6b in the main text, where the phase differences start at x (1) desired at time t = 0, and converge to x desired after a perturbation is applied at t = 50 to force them out from the equilibrium x (1) desired .

A heuristic method to promote stability of functional patterns in positive networks
To promote stability of functional patterns with phase differences |x ij | > π 2 , it is advantageous to design the network weights by minimizing the ones associated to a cos(x ij ) < 0 (i.e., reducing A ij as much as possible) in the cosine-scaled network.Reducing the magnitude of these coupling strengths, or even pruning such interconnections, causes the Gerschgorin disks of the Laplacian L(x desired ) to lie almost entirely in the right half-plane.
In fact, by considering the limit case where all negative connections in the cosine-scaled network are pruned, if the remaining connections describe a connected network, then the Laplacian becomes a classical positive semi-definite Laplacian.This guarantees stability of the desired functional pattern and motivates the following heuristic procedure.
To showcase the effectiveness of the proposed heuristic method to promote stability, we construct an example with n = 7 oscillators, |E| = 9 interconnections, and desired minimum vector of phase differences where x ij = θ j − θ 1 , j = 2, . . ., 7. Notice that the first difference x 12 > π/2, hence cos(x 12 ) < 0. Consider the oscillator network with structural parameters that read as: Such a network admits the following phase-locked equilibrium desired only in x 12 , and has all x ij satisfying |x ij | < π/2 (thus generating a stable functional pattern).Supplementary Fig. 13a-b illustrate the stability of x (0) desired and the functional pattern R 0 associated with this equilibrium.
The network considered in this example has 9 interconnections that can be modified, and the desired equilibrium x (1) desired implies that N = {1}, so that α N is the modification of the first interconnection.To compute the optimal tuning of the network weights we solve the optimization (through standard cvx routines [32]) subject to BD(x)(δ + α) = ω, and (δ + α) ≥ 0.
The optimal α * reads α * = [−0.17060.3152 − 0.8748 59 0.9242 0 0 0 59] T , and the adjusted network adjacency matrix becomes: where the optimal network correction α * has disconnected the entry A 12 (highlighted in red).The matrix Ã guarantees stability of the functional pattern associated with x desired , as we illustrate in Supplementary Fig. 13cd.
In practice, one may jointly optimize the network weights to satisfy the equilibrium constraints and the heuristic stability strategy while trying to keep the overall modification as small as possible.To do so, we let P (resp., N ) denote the set of indices associated with A ij cos(x ij ) > 0 (resp.,A ij cos(x ij ) < 0).Then, the optimization problem that enacts the proposed strategy reads as: where • is a desired vector norm, c 1 , c 2 > 0 are arbitrary penalty coefficients, α 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 tuning vector α that are associated to negative weights δ N in the cosine-scaled network.
To compute the optimal tuning of the network weights in the same 7-oscillator network above, we solve the optimization in Eq. ( 23) with c 1 = 0.1 and c 2 = 10 by minimizing the 1 -norm.The optimal α * reads α * = [−0.17060.3152 − 0.8748 0 0.9242 0 0 0 0] T , and the adjusted network adjacency matrix becomes: , where, as in the procedure above, the optimal network correction α has disconnected the entry A 12 (highlighted in red).Supplementary Fig. 14 illustrates the shift of the Jacobian's eigenvalues while the optimal tuning vector α * is gradually applied.The main differences with respect to the minimization in (21) is that the norm of α * is smaller, and the eigenvalues of Ã are closer to the original ones in A.

Extension of the proposed optimization methods to directed networks
By relaxing the assumption on symmetric adjacency matrices, the matrix form of an oscillator network with The main change in the behavior of directed networks with respect to undirected ones is that the frequencies θ of phase-locked trajectories do not typically converge to the average natural frequency (i.e., θ = ω mean 1).For phase-locked trajectories we have that θ = ω sync 1, where the constant ω sync ∈ R is not known a priori, and can only be estimated in the almost-fully synchronized regime |θ i (t) − θ j (t)| 1 for all i, j ∈ O [69].Yet, not knowing ω sync as we do in the undirected case does not prevent the definition of a framework that can enforce a target functional pattern.In fact, by using equation ( 24) we can concurrently achieve the functional pattern associated with xdesired and assign a desired phase-locked frequency ω sync .To do so, we utilize equation (24) above with θ = ω sync 1 as a constraint in the optimization problems proposed in the main text.This modification allows us to extend any of the control methods to achieve desired functional relationships to directed networks.
As an example of optimization of the coupling strengths, we solve for a network of n = 7 oscillators with adjacency matrix and natural frequencies as follows: , where the entry highlighted in red is the one causing the asymmetry in the network coupling.The target phase differences x 1i are x desired = π T , and the desired phase locking frequency is ω sync = 1 rad/sec.
The solution α * to problem (25) above yields Since the phase time series θ(t) are derived from inherently noisy measurements, we compute the best estimate of the phases θ * (modulo rotation) that are compatible with the noisy measurement in θ(t) by solving the nonconvex phase synchronization problem [11] -that is, the estimation of phases from noisy pairwise relative phase measurements.Given a time window of frequency-synchronized phase time series θ, we find that R ≈ F (see main text and Fig. 16).Moreover, it holds that R − F 2 → 0 as P ij → 1 element-wise.This implies that functional relationships between the phases θ * (t) (encoded in the matrix R) represent the same functional relationships that are measured in fMRI data (encoded in the matrix F ), supporting the usage of Kuramoto oscillators to analyze neural synchronization.

Power network modeling and assumptions
The main manuscript contains an application of our network tuning methods to the IEEE 39 New England power distribution network [7,64].To model the dynamics of this network, we consider a connected power network with generators V 1 and load buses V 2 .A structure-preserving power network model contains |V 1 | second-order Newtonian and |V 2 | first-order kinematic phase oscillators obeying [64]: where M i , D i are the generator inertia constant, and the damping coefficient, respectively.In the equation for the generators V 1 , ω i = P m,i , which is the mechanical power input from the prime mover, and in the equation for the load buses V 2 , ω i = P ,i , which denotes the real power drawn by load i.Finally, the weight a ij equals 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 | is bounded, which is satisfied by frequency-synchronized phase trajectories.Moreover, we assume constant voltage magnitudes |v i | at the loads, so that the weights a ij can be considered fixed.This is a standard assumption in power systems ( also known as decoupling assumption).We refer the interested reader to Ref. [23,Remark 1] for further details.
The generators and bus parameters for the IEEE New England Power network are available in the original article and in classic textbooks [7,64].For simplicity, we set D i = 1 for all loads, which corresponds to a highly damped scenario, possibly due to local excitation controllers.In our simulations, we utilized the standard optimal power flow solver provided by MATPOWER to compute the parameters v, p and θ(0) needed to integrate the Kuramoto model in Matlab through a standard ode45 solver.

Application to additional power network models
Depending on which assumptions are made and on which application is studied, there exist many different power network models in the literature.To demonstrate that the fundamental principles of our procedure remain unchanged even when dealing with more complex dynamics that relax some of the modeling assumptions, we apply our method to two models that differ from the one in equations (26).Henceforth, our goal is to intervene on a power network after a fault occurs between two loads in order to recover a desired (pre-fault) functional pattern.
First, we study the case of a third-order (also known as one-axis) model, such as the one in Ref. [66].The main differences with the structure-preserving power network in equations ( 26) is that the model in Ref. [66] includes the transient dynamics of voltage magnitudes, and that electrical loads are simply modeled as passive impedances.For N = 10 generators, Supplementary Fig. 18a illustrates the reduced power network model obeying the dynamics [66] where θ i is the rotor angle, ω i its frequency, p m,i is the effective mechanical input power of the machine i, M i and D i are the inertia and damping of the mechanical motion, respectively, v i indicates the transient voltage, Im(Y ij ) is the susceptance of the transmission line (i, j), T i denotes the relaxation time of the transient voltage dynamics along the q axis, v f i is the internal voltage, and, finally, χ i and χ i are the static and transient reactances along the d-axis.
Akin to the structure-preserving power network model, stationary operation of the grid corresponds to constant voltages and frequencies, along with constant rotor phase differences.Since our procedure relies on phaselocked trajectories to achieve a desired functional pattern (i.e., power flow), it can be adjusted to intervene on the stationary operation of a challenging model such as the one in equations (27).In fact, in regimes with small voltage and frequency swings, the system is still well approximated by first-order Kuramoto dynamics.Clearly, for a desired functional pattern, the error between the one estimated from the application of our procedure and the one obtained from the dynamics in equations ( 27) will be proportional to the changes in voltages and frequencies.
To test our approach on the model from equations (27), we use the same IEEE 39 New England benchmark case as in the main manuscript.The initial network parameters for our simulations, which represent standard grid operating conditions, are taken from Ref. [77] and Ref. [48].The voltage v i (0) and the initial condition for θ i (0) and ω i (0) = 0 for generator i are fixed using power flow computation.The goal of this application is to recover the initial power flow after the same fault as in Ref. [77] (a line trip between loads 16 and 17) occurs.Such a fault affects the network admittance matrix, and yields an undesired power flow.To recover the pre-fault power flow we follow the same steps as described in the main text (see also Fig. 5b), but because the reduced-network structure is typically a complete graph [21], we do not modify its coupling strengths.Instead, we assume that we can adjust the values of p m ∈ R N , so that p m − diag(D 1 , . . ., D N )ω are the oscillators' natural frequencies that are tuned in order to achieve the pre-fault functional patterns in the post-fault network.The values for P m are computed by rewriting the second one of equations (27) as equation ( 4) in the main text, and solving for the natural frequencies.For the coupling values in δ, we use a ij = |v ss i ||v ss j |Im(Y ij ), where v ss i denotes the steady state value of the voltage v i before applying our intervention.Finally, we sum a positive constant to the obtained natural frequencies so that they are all positive -this constant can be used to adjust the generators to a desired average mechanical input.Supplementary Fig. 18b-c illustrate that our procedure is able to adequately recover the desired functional pattern after the line trips occurs.The error between the pre-fault pattern R 0 and the recovered one R recovered is due to small changes in the frequency and voltage values after P m is modified.We remark that our procedure relies on the dynamics of ω and v leading to small changes.In situations where these changes significantly affect the phases dynamics, the first-order Kuramoto approximation leveraged by our method may not successfully recover a desired functional pattern.
We now turn our attention to models that do not neglect energy losses, such as the one in Ref. [33].By compensating for the losses in the injected power but considering lossy interconnections, the structure-preserving model in equations (26) where ϕ represents the phase shift induced by energy losses.We show in Fig. 19 that our procedure to restore a desired functional pattern in a lossy power network after a fault still works well for small values of the phase shift ϕ.We can recover functional patterns associated to target active power flow conditions by modifying the constraints of our optimization procedures to explicitly take into account the dynamics in equations (28).
At synchronous operating conditions, the dynamics in equations ( 28) reduce to the ones of Kuramoto-Sakaguchi oscillators [63]: The phase shift ϕ makes the matrix formulation of the above dynamics to be incompatible with the ones we have introduced in equation ( 3) of the main text.We can instead write the above dynamics in matrix form by considering the graph G as a directed graph (see also Supplementary Text 1.6 above), with E d being the set of directed edges.That is, we consider each undirected edge as two directed edges where each direction (j, i) = (i, j) has the same weight A ij = A ji .Then, we can write the following equation: where is the diagonal matrix of all sin(x ij + ϕ), and ω sync ∈ R is the synchronization frequency of phase-locked trajectories.
We are now ready to apply the equation ( 29) as a constraint in our numerical optimization routine.Without loss of generality, we set ω sync = 0, and apply the same method developed for the lossless case to the IEEE 39 New England test case to compute the optimal correction after a fault occurs between loads 13 and 14 (the same as in the main text).For a loss ϕ = 0.01, the updated optimization is able to recover the pre-fault functional pattern with a mean error < vec(R) − vec(R recovered ) >= 0.072, where vec(•) vectorizes the matrix (see also Supplementary Fig. 20).This result improves upon the original method proposed for lossless networks by guaranteeing a satisfactory active power flow recovery for losses ϕ that are one order of magnitude larger.
Finally, we observe that for ϕ > 0.01 the fixed parameters of this specific system cause the phases to lose frequency synchronization even at operating conditions.
2 Supplementary Figures    Results of the heuristic method to promote stability of functional patterns containing negative correlations.a Phase differences of the network with adjacency matrix A in Eq. (20).The phase differences converge to x (0) desired .b The functional pattern associated with the phase differences in panel a. c Phase differences of the network with adjacency matrix Ã in Eq. ( 22), after the optimal adjustment is computed from Eq. ( 21).The phase differences converge to x (1) desired .d The functional pattern associated with the phase differences in panel a.Notice that the only rows and columns that change from the functional pattern R 0 are the second row and second column.This is due to the fact that only x 12 = θ 2 − θ 1 differs between the two equilibria x T , where x 12 = θ 2 − θ 1 > π 2 .Notice that, differently from Fig. 7 in the main text, the minimization in equation ( 23) enables a refined optimization of the network weights through the scaling parameters c 1 = 0.1 and c 2 = 10.The left plot illustrates the Gerschgorin disks and the Jacobian's eigenvalues locations for the original network.It can be observed in the zoomed-in panel that one eigenvalue is unstable (λ 2 = 0.0565, in red).The optimal correction α * of the oscillators' coupling strengths is gradually applied from the left-most panel to the right-most one at 1  3 increments.The right zoomed-in panel shows that, as a result of our procedure, n − 1 eigenvalues ultimately lie in the left-hand side of the complex plane (λ 1 = 0 due to rotational symmetry and λ 2 = −0.0178, in green). 2 .The thick black line represents the average 2 -norm distance over 10 4 random initializations, and the shaded area represents the smallest and largest value of the 2 -norm distance.The sudden drop in the largest norm is due to the phases θ ∈ T n evolving in the torus, thus taking values in the interval [0 2π).Our numerical simulations reveal that phase trajectories starting from the other half of the torus (i.e., x 0 − x ∞ > π 2 ) may converge to different stable patterns, which are not compatible with the phases extracted from the functional MRI recordings.b The difference between the functional connectivity F and the functional pattern R. The entries with the largest magnitude are ≈ 0.08, highlighting the stark similarity between the two correlation patterns R and F .The IEEE-39 New England power network reduced to a 10-generator network.Electrical loads are simply modeled as passive impedances.In order to explicitly account for the outside of the system, Generator 1 is assumed to be connected to an infinite bus and has constant phase and frequency [77].b The absolute difference between the pre-fault functional pattern R 0 and the post-fault pattern R fault .The Frobenius norm of this difference is R 0 − R fault F = 0.0907.c The absolute difference between the pre-fault functional pattern R 0 and the recovered pattern R recovered after tuning the power p m at the generators.The Frobenius norm of this difference is R 0 − R recovered F = 0.0653.Error between the pre-fault functional pattern R and the functional pattern R recovered obtained through our procedure as a function of the phase shift ϕ ∈ S 1 .The functional pattern R recovered is computed after a network correction due to a fault that occurs in the IEEE 39 test case between loads 13 and 14 (the same as in the main text).In the presence of a phase shift ϕ, the error between the desired functional pattern and the one associated with the network correction computed by our method remains small for small values of energy loss ϕ.Here, the parameter optimization does not explicitly account for the phase shift ϕ.The mean error is computed as < vec(R) − vec(R recovered ) >, and the maximum error as vec(R) − vec(R recovered ) ∞ , where vec(•) denotes the vectorization.

1 correlationFig 1 .
Fig 1. Network control to enforce a desired functional pattern from an abnormal or undesired one.The left panel contains a network of n = 7 oscillators (top left panel, line thickness is proportional to the coupling strength), whose vector of natural frequencies ω has zero mean.The phase differences with respect to θ 1 (i.e., θ i − θ 1 ) converge to π 8 π 8 π 6 π 6 π 3 2π

34 Fig 2 .
Fig 2.  Mapping between desired phase differences and interconnection weights.a A line network of n = 4 nodes and its parameters.The desired phase differences are shown in red.b Left panel: space of the phase differences; right panel: space of the interconnection weights.The pattern x is illustrated in red in the left panel, and the network weights that achieve such a pattern are represented in red in the right panel.For fixed natural frequencies ω, the green parallelepiped on the left represents the continuum of functional patterns within 0.2 radians from x which can be generated by the positive interconnection weights in the green parallelepiped on the right.

3 A 12 > 0 A 13 > 0 A 23 > 0 A 12 =Fig 3 .
Fig 3.  Algebraic and graph-theoretic conditions for the existence of positive weights that attain a desired functional pattern.a The left side illustrates a simple network of 3 oscillators with adjacency matrix B and vector of natural frequencies ω.The right side illustrates the cone generated by the columns of B. In this example, S = {1, 2} satisfies the conditions for the existence of δ ≥ 0 in equation (5), as ω is contained within the cone generated by the columns B:,S .b The (directed) Hamiltonian path described by the columns of B:,H , with H = {1, 2}, in the network of panel a. c The existence of such an Hamiltonian path, together with a positive projection of ω onto B:,H , also ensure that there exists a strictly positive δ > 0 solution to BD(x)δ = ω.In particular, for any choice of x 12 , x 23 ∈ (0, π), equation (4) reveals that if 0 < A 13 < 0.5/ sin(x 12 + x 23 ), then there exist strictly positive weights A 12 > 0 and A 23 > 0 such that δ > 0.

Fig 7 .
Fig 7. Mechanism underlying the heuristic procedure to promote stability of functional patterns containing negative correlations.For the 7-oscillator network in Supplementary Text 1.5, we apply the procedure in equation (21) to achieve stability of the pattern x desired = 21π 32 π 6 π 6 π 8 π 8 π 3

<Fig 8 .
Fig 8.  Replication of empirically recorded functional connectivity in the brain through tuning of the natural frequencies of Kuramoto oscillators.The anatomical brain organization provides the network backbone over which the oscillators evolve.The filtered fMRI time series provide the relative phase differences between co-fluctuating brain regions, and thus define the desired phase differences x, which is used to calculate the metabolic change encoded in the oscillators' natural frequencies.In this figure, we select the 40-second time window from t 0 = 498 seconds to t f = 538 seconds for subject 18 in the second scanning session.We obtain R − F 2 = 0.2879.Additionally, we verify that the analysis of the Jacobian spectrum (see equation(10)) accurately predicts the stability of the phase-locked trajectories.Supplementary Figure7aillustrates the basin of attraction of R, which we numerically estimate to be half of the torus.

Fig 9 .
Fig 9.  Fault recovery in the IEEE 39 New England power distribution network through minimal and local intervention.a New England power distribution network.The generator terminal buses illustrate the type of generator (coal, nuclear, hydroelectric).We simulate a fault by disconnecting the transmission line 25 (between loads 13 and 14).b The fault causes the voltage phases θ to depart from normal operating conditions, which could 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 and promote a local (sparse) intervention, we solve the optimization in equation (11) by minimizing the 1-norm of the structural parameter modification δ with no scaling parameters in the cost functional.The network returns to the initial operative conditions with a localized modification of the neighboring transmission lines' impedances.

1 † 5 )
thus showing that, in any phase locked trajectory, the oscillators frequency k needs to equal the mean natural frequency 1 n n i=1 ωi.Any feasible functional pattern has n − 1 degrees of freedomThe values of a functional pattern can be uniquely specified using a set of n − 1 correlation values.To see this, let us define the incremental variables x = M θ, where M ∈ R |E|×n is the matrix whose k-th row, associated to xij, is all zeros except for b ki = −1 and b kj = 1.Consider the first n − 1 rows of M , associated to x12, x13, . . .x1n, and notice that they are linearly independent.Moreover, the row associated to xij, i > 1, can be obtained by subtracting the row associated to x1i to the row associated to x1j, implying that the rank of M is n − 1.Any collection of n − 1 linearly independent rows of M defines a full row-rank matrix Mmin (e.g., any n − 1 rows corresponding to the transpose incidence matrix of a spanning tree[30]).We let xmin = Mminθ, where xmin is a smallest set of phase differences that can be used to quantify the synchronization angles among all oscillators.Because ker(Mmin) = 1, we can obtain the phases θ from xmin modulo rotation: θ = M † min xmin − c1, where M † min denotes the Moore-Penrose pseudo-inverse of Mmin and c is some real number.Further, since ker(Mmin) = ker(M ), we can reconstruct all phase differences x from xmin: M M † min xmin = M (θ + c1) = M θ + 0 = x.The above equation reveals that all the differences x are encoded in xmin.That is, any xij can be written as a linear combination of the elements in xmin.For example, if n = 3 and xmin = [x12 x23] T , then x13 is a linear combination of the differences in xmin, i.e., x = xmin, in which x13 = x12 + x23.Thus, because n − 1 incremental variables define all the remaining ones, the entries of any functional pattern must have only n − 1 degrees of freedom.Existence of a strictly positive solution to Problem (Rewrite the pattern assignment problem BD(x)δ = ω as BD(x)δ = B:,HDH,H(x)δH + B :, HD H, H(x)δ H = ω, where the subscripts H and H denote the entries corresponding to the Hamiltonian path in conditions (ii.a)-(ii.b)and the remaining ones, respectively.Since Im( B:,H) = Im(B:,HD:,H(x)) = span(1) ⊥ , Im(B :, HD :, H(x)) ⊆ span(1) ⊥ , and ω ∈ span(1) ⊥ , for any vector δ H, the following set of weights solves the above equation: δH = (B:,HDH,H(x)) † ω − B :, HD H, H(x)δ H = (DH,H(x)B T :,H B:,HDH,H(x)) −1 DH,H(x)B T :,H ω − B :, HD H, H(x)δ H Because the matrix DH,H(x)B T :,H B:,HDH,H(x) is an M-matrix, its inverse has nonnegative entries.Further, by condition (ii.b),DH,H(x)B T in networks with a partition of the oscillators C = {C 1 , . . ., C m }, with C k ⊆ O being a subset of the network oscillators constituting an S 1 -synchronizing graph, k ∈ {1, . . ., m}, where m k=1 C k = O and C k ∩ C = ∅ if k = .Further, ω i = ω j for every i, j ∈ C k , k ∈ {1, . . ., m}. 1 Whenever cluster synchronization emerges, the diagonal blocks (of sizes |C k | × |C k |, k = 1, . . ., m) of the functional pattern associated with partition C satisfy ρ ij = 1 (see Supplementary Fig.

Kuramoto dynamics in equation ( 3 )
of the main text does not hold anymore and requires a rewriting.In what follows, we use the subscript "d" to indicate notation associated to directed graphs.Specifically, E d denote the oriented edge set, so thatD d ∈ R |E d |×|E d | and δ d ∈ R |E d | denote the diagonal matrix of all sin(x ij ) with (j, i) ∈ E d ,and the vector of all the network weights A(i, j) = 0, respectively.Further, let us define B source ∈ R n×|E d | the modified incidence matrix whose columns have nonzero entries only at the edges' sources.That is, B source,k = −1 if k is the source of the interconnection , and B source,k = 0 otherwise.These definitions allow us to define the matrix form for a directed network of Kuramoto oscillators: with v i denoting the nodal voltage magnitude and Y ij being the admittance matrix.The above structure-preserving power network model represents an AC grid with a synchronous generator.Owing to Ref. [23, Lemma 1], the existence and local exponential stability of synchronized solutions of the oscillator model Eq.(26) can be entirely described by means of the first-order Kuramoto model.That is, the load dynamics of a structure-preserving power grid model has the same stable synchronization manifold of Eq. (1) in the main text.

Fig 10 .
Fig 10.Comparison between Pearson correlation coefficient and the metric ρ 12 =< cos(θ 2 − θ 1 ) > t on two phase-locked signals and varying time window lengths.Each point in the plot represents a value of ρ 12 (in blue) and r 12 (in red) computed in a time window [0 T ], where T = 0.1, 0.2, . . ., 5. It can be seen in all panels that ρ 12 is only affected by the phase shift ψ.Instead, the Pearson correlation coefficient returns oscillating values for different window sizes with damping oscillations as the length of the time window increases.a The phase shift in the initial conditions ψ = π 8 .b The phase shift in the initial conditions ψ = π 4 .c The phase shift in the initial conditions ψ = π 3 .d The phase shift in the initial conditions is ψ = π 2 .

Fig 11 .Fig 12 . 6 T
Fig 11.Infinite compatible patterns in a complete graph of identical oscillators.a A network described by a complete graph of n = 6 oscillators with interconnection weights δ = 1 and homogeneous natural frequencies ω = 0. b Each cycle of the network is highlighted in a different color, which is reflected on the entries of the phase differences equilibria x (i) .It holds sin(x (i) ) = [a a a a a a a a b b b b 0 0 0] T ∈ ker(B), where a = sin( π 2 − γ) = sin( π 2 + γ) and b = sin(2γ) = sin(π − 2γ) for all γ ∈ (0, π 2 ).Clearly, any value for γ in the latter interval determines a distinct compatible functional pattern.
Fig 13.Results of the heuristic method to promote stability of functional patterns containing negative correlations.a Phase differences of the network with adjacency matrix A in Eq.(20).The phase differences converge to x

Fig 14 .
Fig 14.Refined mechanism underlying the heuristic procedure to promote stability of functional patterns containing negative correlations.For the 7-oscillator network in Supplementary Text 1.5, we apply the procedure in equation (23) with c 1 = 0.1 and c 2 = 10 to achieve the stability of the pattern x desired = 21π 32 π 6 π 6 π 8 π 8 π 3

4 8 • 10 4 Fig 17 .Fig 18 .
Fig 17.Matrices that describe the network interconnections, the fault, and the local intervention of the power network parameters to recover the pre-fault power distribution.a The adjacency matrix used in the Kuramoto 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 the fault and their immediate neighbors require adjustments.The sparsity of the local intervention is promoted by the usage of the 1 -norm in the optimization problem.

Fig 19 .
Fig 19.Error between the pre-fault functional pattern R and the functional pattern R recovered obtained through our procedure as a function of the phase shift ϕ ∈ S 1 .The functional pattern R recovered is computed after a network correction due to a fault that occurs in the IEEE 39 test case between loads 13 and 14 (the same as in the main text).In the presence of a phase shift ϕ, the error between the desired functional pattern and the one associated with the network correction computed by our method remains small for small values of energy loss ϕ.Here, the parameter optimization does not explicitly account for the phase shift ϕ.The mean error is computed as < vec(R) − vec(R recovered ) >, and the maximum error as vec(R) − vec(R recovered ) ∞ , where vec(•) denotes the vectorization.

1 Fig 20 .
Fig 20.Nominal, post-fault, and recovered functional pattern in a network with lossy communications.In all panels, the loss is fixed to ϕ = 0.01.a Functional pattern associated to nominal power flow in the IEEE 39 New England test case.b Functional pattern associated to a power flow disruption due to a fault that disconnects loads 13 and 14. c The recovered functional pattern after our procedure with updated constraint (equation (29)) is applied.The mean error between the pre-fault and the recovered functional patterns is < vec(R) − vec(R recovered ) >= 0.072, where vec(•) denotes the vectorization of the patterns.
the active power load at node i, where D i is the damping coefficient, andA ij = |v i ||v j |I(Y ij )/D i ,with v i denoting the nodal voltage magnitude and I(Y ij ) being the imaginary part of the admittance Y ij (see Supplementary Information for details about this model).In this example, we choose a highly damped scenario