# Functional Control of Oscillator Networks

Tommaso Menara  Department of Mechanical and Aerospace Engineering, University
of California, San Diego, La Jolla, CA 92093, USA  Giacomo Baggio  Department
of Information Engineering, University of Padova, Padova, 35131, Italy
Danielle S. Bassett  Department of Physics & Astronomy, Department of
Bioengineering, Department of Electrical & Systems Engineering, Department of
Neurology, Department of Psychiatry, University of Pennsylvania, Philadelphia,
PA 19104, USA
The Santa Fe Institute, Santa Fe, NM 87506, USA  Fabio Pasqualetti  Department
of Mechanical Engineering, University of California, Riverside, Riverside, CA
92521, USA<EMAIL_ADDRESS>
## Abstract

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 [17, 58, 15]. For
instance, distinctive network-wide patterns of synchrony [72, 31, 71]
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
$\mathcal{G}=\\{\mathcal{O},\mathcal{E}\\}$ of $n$ oscillators with dynamics

$\displaystyle\dot{\theta}_{i}=\omega_{i}+\sum_{j=1}^{n}A_{ij}\sin(\theta_{j}-\theta_{i}),$
(1)

where $\omega_{i}\in\mathbb{R}$ and $\theta_{i}\in\mathbb{S}^{1}$ are the
frequency and phase of the $i$-th oscillator, respectively, $A=[A_{ij}]$ is
the weighted adjacency matrix of $\mathcal{G}$, and
$\mathcal{O}=\\{1,\dots,n\\}$ and
$\mathcal{E}\subseteq\mathcal{O}\times\mathcal{O}$ denote the oscillator and
interconnection sets, respectively. In this work we consider the case where
the network $\mathcal{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 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 $\theta_{i}(t)$ and $\theta_{j}(t)$, we define the
correlation coefficient

$\rho_{ij}=<\cos(\theta_{j}(t)-\theta_{i}(t))>_{t},$ (2)

where $<\cdot>_{t}$ denotes the average over time. A _functional pattern_ is
formally defined as the symmetric matrix $R$ whose $i,j$-th entry equals
$\rho_{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, $\rho_{ij}\to 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.

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 ${\boldsymbol{\omega}}$ has
zero mean. The phase differences with respect to $\theta_{1}$ (i.e.,
$\theta_{i}-\theta_{1}$) converge to
$\left[\frac{\pi}{8}~{}\frac{\pi}{8}~{}\frac{\pi}{6}~{}\frac{\pi}{6}~{}\frac{\pi}{3}~{}\frac{2\pi}{3}\right]$,
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 $\boldsymbol{\omega}$ are adjusted to $A_{\mathrm{c}}$ and
${\boldsymbol{\omega}}_{\mathrm{c}}$) to obtain the phase differences
$\left[\frac{2\pi}{3}~{}\frac{\pi}{3}~{}\frac{\pi}{6}~{}\frac{\pi}{6}~{}\frac{\pi}{8}~{}\frac{\pi}{8}\right]$,
which encode the desired functional pattern in the bottom right. In this
example, we have computed the closest set (in the $\ell_{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.

While synchronization phenomena in oscillator networks have been studied
extensively (e.g., see [57, 74, 70, 52, 10]), the development of control
methods to impose desired synchronous behaviors has only recently attracted
the attention of the research community [55, 25, 28, 45]. 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.

## Analytical results

### Feasible functional patterns in positive networks

A functional pattern is an $n\times 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, 81, 39].

For the undirected network $\mathcal{G}=(\mathcal{O},\mathcal{E})$, let
$x_{ij}=\theta_{j}-\theta_{i}$ be the difference of the phases of the
oscillators $i$ and $j$, and let $\boldsymbol{x}\in\mathbb{R}^{|\mathcal{E}|}$
be the vector of all phase differences with $(i,j)\in\mathcal{E}$ and $i<j$.
The network dynamics (1) can be conveniently rewritten in vector form as (see
Methods)

$\displaystyle
BD({\boldsymbol{x}})\boldsymbol{\delta}=[\omega_{1}\;~{}\cdots\;~{}\omega_{n}]^{\mathsf{T}}-\dot{\boldsymbol{\theta}},$
(3)

where $B\in\mathbb{R}^{n\times|\mathcal{E}|}$ is the (oriented) incidence
matrix of the network $\mathcal{G}$,
$D({\boldsymbol{x}})\in\mathbb{R}^{|\mathcal{E}|\times|\mathcal{E}|}$ is a
diagonal matrix of the sine functions in equation (1), and
$\boldsymbol{\delta}\in\mathbb{R}^{|\mathcal{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
$\dot{\boldsymbol{\theta}}$ satisfies
$\dot{\boldsymbol{\theta}}=\omega_{\text{mean}}\textbf{1}$, where
$\omega_{\text{mean}}=\left(\frac{1}{n}\sum_{i=1}^{n}\omega_{i}\right)$ is the
average of the natural frequencies of the oscillators. Further, since
$\mathcal{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},\dots,x_{n-1,n}\\}$. In fact, for any pair of
oscillators $i$ and $j$ with $i<j$, it holds
$x_{ij}=\sum_{k=i}^{j-1}x_{k,k+1}$. This implies that the vector of all phase
differences in equation (3), and in fact any $n\times n$ functional pattern,
has only $n-1$ degrees of freedom and can be uniquely specified with a set of
$n-1$ independent differences $\boldsymbol{x}_{\text{desired}}$ (see Methods).
Following this reasoning and to avoid cluttered notation, let
$\boldsymbol{\omega}=[\omega_{1}-\omega_{\text{mean}}\;~{}\cdots\;~{}\omega_{n}-\omega_{\text{mean}}]^{\mathsf{T}}$,
and notice that the problem of enforcing a desired functional pattern
simplifies to (i) converting the desired functional pattern to the
corresponding phase differences $\boldsymbol{x}_{\text{desired}}$, and (ii)
computing the network weights $\boldsymbol{\delta}$ to satisfy the following
equation:

$BD({\boldsymbol{x}})\boldsymbol{\delta}=\boldsymbol{\omega},$ (4)

where we note that the vector $\boldsymbol{\omega}$ has zero mean and that,
with a slight abuse of notation, $D({\boldsymbol{x}})$ denotes the
$|\mathcal{E}|$-dimensional diagonal matrix of the sine of the phase
differences uniquely defined by the ($n-1$)-dimensional vector
$\boldsymbol{x}_{\text{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 $\boldsymbol{x}$,
this problem reads as

$\displaystyle\mathrm{find}$ $\displaystyle\boldsymbol{\delta}$ (5)
$\displaystyle\mathrm{subject~{}to}$ $\displaystyle
BD(\boldsymbol{x})\boldsymbol{\delta}=\boldsymbol{\omega},$ (5a) and
$\displaystyle\boldsymbol{\delta}\geq 0.$ (5b)

It should be noticed that the feasibility of the optimization problem (5)
depends on the sign of the entries of the diagonal matrix $D(\boldsymbol{x})$,
but is independent of their magnitude. To see this, notice that

$\displaystyle
D(\boldsymbol{x})=\text{sign}(D(\boldsymbol{x}))|D(\boldsymbol{x})|,$

where the $\mathrm{sign}(\cdot)$ and absolute value $|\cdot|$ operators are
applied element-wise. Then, Problem (5) is feasible if and only if there
exists a nonnegative solution to

$\displaystyle\underbrace{B\text{sign}(D(\boldsymbol{x}))}_{\bar{B}}\underbrace{|D(\boldsymbol{x})|\boldsymbol{\delta}}_{\bar{\boldsymbol{\delta}}}=\boldsymbol{\omega}.$

The feasibility of the latter equation, in turn, depends on the projections of
the natural frequencies $\boldsymbol{\omega}$ on the columns of $\bar{B}$: a
nonnegative solution exists if $\boldsymbol{\omega}$ belongs to the cone
generated by the columns of $\bar{B}$. This also implies that, if a network
admits a desired functional pattern $\boldsymbol{x}$ then, by tuning its
weights, the same network can generate any other functional pattern
${\boldsymbol{x}}_{\text{new}}$ such that
$\text{sign}(D({\boldsymbol{x}}_{\text{new}}))=\text{sign}(D({\boldsymbol{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.

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
${\boldsymbol{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 $\boldsymbol{\omega}$, the green parallelepiped
on the left represents the continuum of functional patterns within $0.2$
radians from ${\boldsymbol{x}}$ which can be generated by the positive
interconnection weights in the green parallelepiped on the right.

A sufficient condition for the feasibility of Problem (5) is as follows:

_There exists $\boldsymbol{\delta}\geq 0$ such that
$BD(\boldsymbol{x})\boldsymbol{\delta}={\boldsymbol{\omega}}$ if there exists
a set $\mathcal{S}$ satisfying:_

  1. (i.a)

_$D_{ii}(\boldsymbol{x})D_{jj}(\boldsymbol{x})B_{:,i}^{\mathsf{T}}B_{:,j}\leq
0$ for all $i,j\in\mathcal{S}$ with $i\neq j$ and $D_{ii},D_{jj}\neq 0$;_

  2. (i.b)

_${\boldsymbol{\omega}}^{\mathsf{T}}B_{:,i}D_{ii}(\boldsymbol{x}) >0$ for all
$i\in\mathcal{S}$;_

  3. (i.c)

_${\boldsymbol{\omega}}\in\mathrm{Im}(B_{:,\mathcal{S}})$._

Equivalently, the above conditions ensure that ${\boldsymbol{\omega}}$ is
contained within the cone generated by the columns of
$\bar{B}_{:,\mathcal{S}}$ (see Fig. 3a for a self-contained example). To see
this, rewrite the pattern assignment problem
$BD({\boldsymbol{x}})\boldsymbol{\delta}={\boldsymbol{\omega}}$ as

$\displaystyle
BD({\boldsymbol{x}})\boldsymbol{\delta}=B_{:,\mathcal{S}}D_{\mathcal{S},\mathcal{S}}({\boldsymbol{x}})\boldsymbol{\delta}_{\mathcal{S}}+B_{:,\tilde{\mathcal{S}}}D_{\tilde{\mathcal{S}},\tilde{\mathcal{S}}}({\boldsymbol{x}})\boldsymbol{\delta}_{\tilde{\mathcal{S}}}=\boldsymbol{\omega},$
(6)

where the subscripts $\mathcal{S}$ and $\tilde{\mathcal{S}}$ denote the
entries corresponding to the set $\mathcal{S}$ and the remaining ones,
respectively. If the vectors $B_{:,i}$, $i\in{\mathcal{S}}$, are linearly
independent, condition (i.a) implies that
$D_{\mathcal{S},\mathcal{S}}B_{:,\mathcal{S}}^{\mathsf{T}}B_{:,{\mathcal{S}}}D_{\mathcal{S},\mathcal{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 $\mathcal{S}$ with any subset
$\mathcal{S}_{\text{ind}}\subset\mathcal{S}$ such that the vectors $B_{:,i}$,
$i\in{\mathcal{S}_{\text{ind}}}$, are linearly independent. Condition (i.c)
guarantees the existence of a solution to
$B_{:,\mathcal{S}}D_{\mathcal{S},\mathcal{S}}\boldsymbol{\delta}_{\mathcal{S}}=\boldsymbol{\omega}$.
A particular solution to the latter equation is

$\boldsymbol{\delta}_{\mathcal{S}}=\left(B_{:,\mathcal{S}}D_{\mathcal{S},\mathcal{S}}({\boldsymbol{x}})\right)^{\dagger}\boldsymbol{\omega}=(D_{\mathcal{S},\mathcal{S}}({\boldsymbol{x}})B_{:,\mathcal{S}}^{\mathsf{T}}B_{:,\mathcal{S}}D_{\mathcal{S},\mathcal{S}}({\boldsymbol{x}}))^{-1}D_{\mathcal{S},\mathcal{S}}({\boldsymbol{x}})B_{\mathcal{S},\mathcal{S}}^{\mathsf{T}}\boldsymbol{\omega}>0$

where $(\cdot)^{\dagger}$ denotes the Moore-Penrose pseudo-inverse of a
matrix. The positivity of $\boldsymbol{\delta}_{\mathcal{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
$\boldsymbol{\delta}_{\mathcal{S}}=\left(B_{:,\mathcal{S}}D_{\mathcal{S},\mathcal{S}}({\boldsymbol{x}})\right)^{\dagger}\boldsymbol{\omega}>0$
and $\boldsymbol{\delta}_{\tilde{\mathcal{S}}}=\boldsymbol{0}$.

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 $\bar{B}$ and vector
of natural frequencies ${\boldsymbol{\omega}}$. The right side illustrates the
cone generated by the columns of $\bar{B}$. In this example,
$\mathcal{S}=\\{1,2\\}$ satisfies the conditions for the existence of
$\boldsymbol{\delta}\geq 0$ in equation (5), as ${\boldsymbol{\omega}}$ is
contained within the cone generated by the columns $\bar{B}_{:,\mathcal{S}}$.
b The (directed) Hamiltonian path described by the columns of
$\bar{B}_{:,\mathcal{H}}$, with $\mathcal{H}=\\{1,2\\}$, in the network of
panel a. c The existence of such an Hamiltonian path, together with a positive
projection of ${\boldsymbol{\omega}}$ onto $\bar{B}_{:,\mathcal{H}}$, also
ensure that there exists a strictly positive $\boldsymbol{\delta}>0$ solution
to $BD({\boldsymbol{x}})\boldsymbol{\delta}={\boldsymbol{\omega}}$. In
particular, for any choice of $x_{12},x_{23}\in{(0,~{}\pi)}$, 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 $\boldsymbol{\delta}>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, $\boldsymbol{\delta}>0$ rather than
$\boldsymbol{\delta}\geq 0$). While, in general, this is a considerably harder
problem, a sufficient condition for the existence of a strictly positive
solution $\boldsymbol{\delta}>0$ is that the network with incidence matrix
$\bar{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 $\boldsymbol{\delta}>0$ to
$BD({\boldsymbol{x}}){\boldsymbol{\delta}}=\boldsymbol{\omega}$ if_

  * (ii.a)

_the network with incidence matrix $\bar{B}$ contains a directed Hamiltonian
path $\mathcal{H}$;_

  * (ii.b)

_${\boldsymbol{\omega}}^{\mathsf{T}}B_{:,i}D_{ii}(\boldsymbol{x}) >0$ for all
$i\in\mathcal{H}$;_

The incidence matrix $\bar{B}_{:,\mathcal{H}}$ of a directed Hamiltonian path
$\mathcal{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
$\boldsymbol{\omega}\in\mathrm{Im}(B_{:,\mathcal{H}})$. Second, the columns of
the incidence matrix $\bar{B}_{:,\mathcal{H}}$ satisfy
$\bar{B}_{:,i}^{\mathsf{T}}\bar{B}_{:,i}=2$ and
$\bar{B}_{:,i}^{\mathsf{T}}\bar{B}_{:,j}\in\\{0,-1\\}$ for all
$i,j\in\mathcal{H}$, $i\neq j$. Then, letting the set $\mathcal{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 $\boldsymbol{\delta}$ that solves the pattern assignment problem
$BD({\boldsymbol{x}})\boldsymbol{\delta}={\boldsymbol{\omega}}$. Furthermore,
by rewriting the pattern assignment problem as in equation (6), the following
vector of interconnection weights solves such equation (see Methods):

$\displaystyle\boldsymbol{\delta}_{\mathcal{H}}$
$\displaystyle=\left(B_{:,\mathcal{H}}D_{\mathcal{H},\mathcal{H}}({\boldsymbol{x}})\right)^{\dagger}\left(\boldsymbol{\omega}-B_{:,\tilde{\mathcal{H}}}D_{\tilde{\mathcal{H}},\tilde{\mathcal{H}}}({\boldsymbol{x}})\boldsymbol{\delta}_{\tilde{\mathcal{H}}}\right).$

Because $\bar{B}_{:,\mathcal{H}}=B_{:,\mathcal{H}}D_{\mathcal{H},\mathcal{H}}$
defines an Hamiltonian path and because of (ii.b), the vector
$\left(B_{:,\mathcal{H}}D_{\mathcal{H},\mathcal{H}}({\boldsymbol{x}})\right)^{\dagger}\boldsymbol{\omega}$
contains only strictly positive entries. Thus, for any sufficiently small and
positive vector $\boldsymbol{\delta}_{\tilde{\mathcal{H}}}$, the weights
$\boldsymbol{\delta}_{{\mathcal{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 $\omega_{i}<\omega_{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 $\frac{\pi}{2}$ (trivially, when the phase differences
are smaller than $\frac{\pi}{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 $\boldsymbol{\delta}$, 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., $\boldsymbol{\delta}$ with nonzero components) all
compatible equilibria ${\boldsymbol{x}}^{(i)}$, $i\in\\{1,\dots,\ell\\}$, must
satisfy

$D(\boldsymbol{x}^{(i)})\boldsymbol{\delta}=B^{\dagger}\boldsymbol{\omega}+\mathrm{ker}(B).$
(7)

From equation (7), we can see that the sine vector of all compatible
equilibria must belong to a specific affine subspace of
$\mathbb{R}^{|\mathcal{E}|}$:

$\sin(\boldsymbol{x}^{(i)})\in\mathrm{diag}(\boldsymbol{\delta})^{-1}\left(B^{\dagger}\boldsymbol{\omega}+\mathrm{ker}(B)\right).$
(8)

Rewriting equation (4) in the above form connects the existence of distinct
functional patterns with $\mathrm{ker}(B)$, the latter featuring a number of
well known properties. For instance, it holds that
$\mathrm{dim}(\mathrm{ker}(B))=|\mathcal{E}|-n+c$, where $c$ is the number of
connected components in a network, and that $\mathrm{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 $\boldsymbol{x}$, we can write
$\sin({\boldsymbol{x}})=[\sin(\boldsymbol{x}_{\mathrm{desired}}^{\mathsf{T}})~{}\sin(\boldsymbol{x}_{\mathrm{dep}}^{\mathsf{T}})]^{\mathsf{T}}$,
where $\boldsymbol{x}_{\mathrm{dep}}$ denotes the phase differences dependent
on $n-1$ desired phase differences $\boldsymbol{x}_{\mathrm{desired}}$. Thus,
all the $\boldsymbol{x}_{\mathrm{desired}}$ for which
$\sin(\boldsymbol{x}_{\mathrm{dep}})$ intersects the affine space described by
$\mathrm{diag}(\boldsymbol{\delta})^{-1}\left(B^{\dagger}\boldsymbol{\omega}+\mathrm{ker}(B)\right)$
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 $\boldsymbol{\delta}=\boldsymbol{1}$ and
$\boldsymbol{\omega}=\boldsymbol{0}$, so that equation (8) holds whenever
$\sin(\boldsymbol{x}^{(i)})\in\mathrm{ker}(B)$. In networks with tree
topologies it holds $\mathrm{ker}(B)=\boldsymbol{0}$, and
$\sin(\boldsymbol{x}^{(i)})=\boldsymbol{0}$ is satisfied by $2^{n-1}$ patterns
of the form $x^{(i)}_{jk}=0,\pi$, for all $(j,k)\in\mathcal{E}$. Consider now
cycle networks, where $\mathrm{ker}(B)=\mathrm{span}~{}\boldsymbol{1}$. For
any cycle of $n\geq 3$ oscillators, two families of patterns are
straightforward to derive. First, there are $2^{n-1}$ patterns of the form
$\boldsymbol{x}^{(i)}_{k,k+1}=0,\pi$, with $k=1,\dots,n-1$, and
$x^{(i)}_{n1}=-\sum_{k=1}^{n-1}x^{(i)}_{k,k+1}$. Second, there are $n-1$ splay
states [22], where the oscillators’ phases evenly span the unitary circle,
with $x^{(i)}_{jk}=\frac{2\pi m}{n}$, $m=1,\dots,n-1$, $(j,k)\in\mathcal{E}$.
Fig. 4 illustrates the compatible functional patterns satisfying equation (8)
in a positive network of three fully synchronizing oscillators.

Fig 4:  The intersection of an affine space with
$\sin(\boldsymbol{x}_{\mathrm{dep}})$ determines the compatible functional
patterns of 3 identical oscillators. Consider a fully connected network of
$n=3$ identical oscillators with zero natural frequency and
$\boldsymbol{\delta}=\boldsymbol{1}$. It is well known that
${\boldsymbol{x}}^{(1)}=[0~{}0]^{\mathsf{T}}$,
${\boldsymbol{x}}^{(2)}=[\pi~{}0]^{\mathsf{T}}$,
${\boldsymbol{x}}^{(3)}=[0~{}\pi]^{\mathsf{T}}$,
${\boldsymbol{x}}^{(4)}=[\pi~{}\pi]^{\mathsf{T}}$ are phase difference
equilibria. Furthermore, because $\sin(\theta)=\sin(\pi-\theta)$, this figure
illustrates $\sin(x_{13})$ as a function of $x_{12}$ and $x_{23}$ in four
different panels: $\sin(x_{12}+x_{23})$ (top left), $\sin(\pi-x_{12}+x_{23})$
(top right), $\sin(x_{12}+\pi-x_{23})$ (bottom left), and
$\sin(-x_{12}-x_{23})$ (bottom right). The fourth panel reveals that the two
functional patterns compatible with ${\boldsymbol{x}}^{(j)}$,
$j\in\\{1,\dots,4\\}$, correspond to
${\boldsymbol{x}}^{(5)}=[2\pi/3~{}2\pi/3]^{\mathsf{T}}$ and
${\boldsymbol{x}}^{(6)}=[-2\pi/3~{}-2\pi/3]^{\mathsf{T}}$ (in red).

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 $\gamma\in\mathbb{S}^{1}$ in a cycle of
$n=4$ oscillators. Finally, as complete graphs 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 $\ell\geq 1$ desired functional patterns
$\\{{\boldsymbol{x}}^{(i)}\\}_{i=1}^{\ell}$. We first notice that equation (7)
reveals that to achieve a desired functional pattern ${\boldsymbol{x}}^{(i)}$
with components not equal to $k\pi$, $k\in\mathbb{Z}$, the network weights
$\boldsymbol{\delta}$ must belong to an $|\mathcal{E}|$-dimensional affine
subspace of $\mathbb{R}^{|\mathcal{E}|}$:

$\boldsymbol{\delta}\in
D(\boldsymbol{x}^{(i)})^{-1}\left(B^{\dagger}\boldsymbol{\omega}+\mathrm{ker}(B)\right),\quad\forall
i=1,\dots,\ell.$ (9)

Let
$\Gamma_{i}=D(\boldsymbol{x}^{(i)})^{-1}\left(B^{\dagger}\boldsymbol{\omega}+\mathrm{ker}(B)\right)$.
Then, to concurrently realize a collection of patterns
$\\{{\boldsymbol{x}}^{(i)}\\}_{i=1}^{\ell}$, a solution to equations (9)
exists if and only if $\bigcap_{i=1}^{\ell}\Gamma_{i}\neq\emptyset$. It is
worth noting that, whenever the latter intersection coincides with a
singleton, then there exists a single choice of network weights that realizes
$\\{{\boldsymbol{x}}^{(i)}\\}_{i=1}^{\ell}$. However, if
$\bigcap_{i=1}^{\ell}\Gamma_{i}$ corresponds to a subspace, then infinite
networks can realize the desired collection of functional patterns. We
conclude by emphasizing that a positive $\boldsymbol{\delta}$ that achieves
the desired patterns exists if and only if
$\left(\bigcap_{i=1}^{\ell}\Gamma_{i}\right)\cap\mathbb{R}_{\geq
0}^{|\mathcal{E}|}\neq\emptyset$. That is, if the network weights belong to
the nonempty intersection of the $\ell$ affine subspaces with the positive
orthant.

Fig 5:  A homogeneous cycle network admits infinite compatible functional
patterns. Since $\mathrm{ker}(B)=\mathrm{span}~{}\boldsymbol{1}$, the cycle
network admits infinite compatible equilibria, which can be parameterized by
$\gamma\in\mathbb{S}^{1}$ as
$\boldsymbol{x}^{(i)}(\gamma)=[\frac{\pi}{2}-\gamma,~{}\frac{\pi}{2}+\gamma,~{}\frac{\pi}{2}-\gamma,~{}\frac{\pi}{2}+\gamma]^{\mathsf{T}}$.
Any arbitrarily small variation of $\gamma$ yields
$\sin(\boldsymbol{x}^{(i)}(\gamma))\in\mathrm{ker}(B)$. The right panel
illustrates the patterns associated with $\boldsymbol{x}^{(i)}(\gamma)$,
$i=1,\dots,5$ for increments of $\gamma$ of 0.2 radians.

### 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]

$J=\frac{\partial}{\partial\boldsymbol{\theta}}\dot{\boldsymbol{\theta}}=-\underbrace{B\operatorname{diag}\left(\\{A_{ij}\cos({x_{ij}})\\}_{(i,j)\in\mathcal{E}}\right)B^{\mathsf{T}}}_{\mathcal{L}({\boldsymbol{x}})},$
(10)

where $\mathcal{L}({\boldsymbol{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(\theta_{j}-\theta_{i})$). The
functional pattern ${\boldsymbol{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 $\frac{\pi}{2}$ (that is, the infinity-
norm of $\boldsymbol{x}$ satisfies
$\|{\boldsymbol{x}}\|_{\infty}<\frac{\pi}{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 $\boldsymbol{\delta}$ such that
$A_{ij}>0$ if $|x_{ij}|<\frac{\pi}{2}$ and $A_{ij}<0$ otherwise, so that the
matrix $\mathcal{L}$ becomes the Laplacian of a positive network (see
Methods). Furthermore, we observe that in the particular case where some
differences $|x_{ij}|=\frac{\pi}{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 $\frac{\pi}{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
$\mathcal{L}_{ij}>0$ whenever $|x_{ij}|>\frac{\pi}{2}$, thus possibly changing
the sign of its diagonal entries $\mathcal{L}_{ii}=\sum_{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 $\mathcal{L}$ is structurally balanced if and only if its oscillators
can be partitioned into two sets, $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$,
such that all $(i,j)\in\mathcal{E}$ with $A_{ij}\cos(x_{ij})<0$ connect
oscillators in $\mathcal{O}_{1}$ to oscillators in $\mathcal{O}_{2}$, and all
$(i,j)\in\mathcal{E}$ with $A_{ij}\cos(x_{ij})>0$ connect oscillators within
$\mathcal{O}_{i}$, $i\in\\{1,2\\}$. If a network is structurally balanced,
then its Laplacian has mixed eigenvalues [86]. Therefore, we conclude the
following:

_If the functional pattern ${\boldsymbol{x}}$ yields a structurally balanced
cosine-scaled network, then ${\boldsymbol{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, $\boldsymbol{x}$ is unstable
whenever $|x_{ij}|>\frac{\pi}{2}$ for any $i,j$. Instead, for a cycle network
with positive weights, the pattern $\boldsymbol{x}$ can be stable only if it
contains at most one phase difference $\frac{\pi}{2}<|x_{ij}|<\gamma$, where
$\gamma\approx 1.789776$ solves $\gamma-\tan(\gamma)=2\pi$ (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 ${\boldsymbol{x}}$ and network weights $\boldsymbol{\delta}$, we seek
to solve

$\displaystyle\min_{\boldsymbol{\alpha}}$
$\displaystyle\left\|\boldsymbol{\alpha}\right\|_{2}$ (11)
$\displaystyle\mathrm{subject~{}to}$ $\displaystyle
BD({\boldsymbol{x}})(\boldsymbol{\delta}+\boldsymbol{\alpha})={\boldsymbol{\omega}},$
(11a) and $\displaystyle(\boldsymbol{\delta}+\boldsymbol{\alpha})\geq 0,$
(11b)

where $\boldsymbol{\alpha}\in\mathbb{R}^{|\mathcal{E}|}$ are the controllable
modifications of the network weights, and $\|\cdot\|_{2}$ denotes the
$\ell^{2}$-norm. Fig. 6a illustrates the control of a functional pattern in a
line network of $n=4$ oscillators.

Fig 6: Optimal interventions for desired functional patterns. a For the line
network in Fig. 2a, we solve Problem (11) to assign the desired pattern
$\boldsymbol{x}_{\text{desired}}=[\frac{4\pi}{5}~{}\frac{\pi}{3}~{}\frac{\pi}{10}]^{\mathsf{T}}$.
The starting pattern
$\boldsymbol{x}_{\text{original}}=[\frac{\pi}{10}~{}\frac{\pi}{3}~{}\frac{4\pi}{5}]^{\mathsf{T}}$
is associated to interconnection weights
$\boldsymbol{\delta}=[3.4026~{}3.4641~{}6.4721]^{\mathsf{T}}$. Applying the
optimal correction $\boldsymbol{\alpha}^{*}$ yields positive interconnection
weights
$\boldsymbol{\delta}+\boldsymbol{\alpha}^{*}=[6.4721~{}3.4026~{}3.4641]^{\mathsf{T}}$
that achieve the desired functional patterns
$\boldsymbol{x}_{\text{desired}}$. b Joint allocation of two compatible
equilibria for the phase difference dynamics. By taking $\theta_{1}$ as a
reference, we choose two points for the phase differences
$x_{1i}=\theta_{i}-\theta_{1}$, $i\in\\{2,\dots,7\\}$, to be imposed as
equilibria in a network of $n=7$ oscillators:
${\boldsymbol{x}}^{(1)}_{\text{desired}}=\left[\frac{\pi}{6}~{}\frac{\pi}{6}~{}\frac{\pi}{4}~{}\frac{\pi}{4}~{}\frac{\pi}{6}~{}\frac{\pi}{4}\right]^{\mathsf{T}}$
and
${\boldsymbol{x}}^{(2)}_{\text{desired}}=\left[\frac{\pi}{8}~{}\frac{\pi}{3}~{}\frac{\pi}{4}~{}\frac{\pi}{4}~{}\frac{\pi}{6}~{}\frac{\pi}{4}\right]^{\mathsf{T}}$.
In this example, we find a set of interconnection weights
$(\boldsymbol{\delta}+\boldsymbol{\alpha}^{*})$ that solves the minimization
problem (11) with constraint (12). The trajectories start at the (unstable)
equilibrium point ${\boldsymbol{x}}^{(1)}_{\text{desired}}$ at time $t=0$, and
converge to the point ${\boldsymbol{x}}^{(2)}_{\text{desired}}$ after a small
perturbation $\boldsymbol{p}\in\mathbb{T}^{7}$, with $p_{i}\in[0~{}0.05]$, is
applied to the phase difference trajectories at time $t=50$.

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 $\\{\boldsymbol{x}^{(i)}\\}_{i=1}^{\ell}$. To do so, we
simply replace the constraint (11a) with

$BD(\boldsymbol{x}^{(i)})(\boldsymbol{\delta}+\boldsymbol{\alpha})=\boldsymbol{\omega},\quad\forall
i=1,\dots,\ell.$ (12)

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 $\boldsymbol{x}$ is stable for the network with weights
$\boldsymbol{\delta}+\boldsymbol{\alpha}^{*}$. To promote stability of the
pattern $\boldsymbol{x}$, we use a heuristic procedure based on the classic
Gerschgorin’s theorem [29]. Recall that stability of $\boldsymbol{x}$ is
guaranteed when the associated Jacobian matrix has a Laplacian structure, with
negative diagonal entries and nonnegative off-diagonal entries. Further,
instability of $\boldsymbol{x}$ depends primarily on the negative off-diagonal
entries $A_{ij}\cos(x_{ij})$ of the Jacobian (these entries 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
$\boldsymbol{\delta}_{\mathcal{N}}$ and $\boldsymbol{\alpha}_{\mathcal{N}}$
denote the entries of the weights $\boldsymbol{\delta}$ and tuning vector
$\boldsymbol{\alpha}$, 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:

$\displaystyle\min_{\boldsymbol{\alpha}}$ $\displaystyle\ \
\left\|\boldsymbol{\delta}_{\mathcal{N}}+\boldsymbol{\alpha}_{\mathcal{N}}\right\|_{2}$
(13) $\displaystyle\mathrm{subject~{}to}$ $\displaystyle\ \
BD({\boldsymbol{x}})(\boldsymbol{\delta}+\boldsymbol{\alpha})={\boldsymbol{\omega}},$
$\displaystyle\mathrm{and}$ $\displaystyle\ \
(\boldsymbol{\delta}+\boldsymbol{\alpha})\geq 0,$

Carefully reducing the weights
$\boldsymbol{\delta}_{\mathcal{N}}+\boldsymbol{\alpha}_{\mathcal{N}}$ promotes
stability of the target pattern. Fig. 7 illustrates the shift of the
Jacobian’s eigenvalues while the optimal tuning vector
$\boldsymbol{\alpha}^{*}$ is gradually applied to a $7$-oscillator network to
achieve stability of a functional pattern containing negative correlations
(the network parameters can be found in the Supplementary Information).
Finally, we remark that the procedure in (21) can be further refined by
introducing scaling constants to penalize
$\left\|\boldsymbol{\delta}_{\mathcal{N}}+\boldsymbol{\alpha}_{\mathcal{N}}\right\|_{2}$
differently from the modification of other interconnection weights (see
Supplementary Information for further details and an example).

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
${\boldsymbol{x}}_{\text{desired}}=\left[\frac{21\pi}{32}~{}\frac{\pi}{6}~{}\frac{\pi}{6}~{}\frac{\pi}{8}~{}\frac{\pi}{8}~{}\frac{\pi}{3}\right]^{\mathsf{T}}$,
where $x_{12}=\theta_{2}-\theta_{1}>\frac{\pi}{2}$. 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 ($\lambda_{2}=0.0565$, in red). The optimal correction
$\boldsymbol{\alpha}^{*}$ is gradually applied to the existing
interconnections from the left-most panel to the right-most one at
$\frac{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 ($\lambda_{1}=0$ due to rotational symmetry and
$\lambda_{2}=-0.0178$, in green).

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

$\displaystyle\min_{\boldsymbol{\alpha},\boldsymbol{\beta}}$
$\displaystyle\left\|\begin{bmatrix}\boldsymbol{\alpha}^{\mathsf{T}}&\boldsymbol{\beta}^{\mathsf{T}}\end{bmatrix}\right\|_{2}$
(14) $\displaystyle\mathrm{subject~{}to}$ $\displaystyle
BD({\boldsymbol{x}})(\boldsymbol{\delta}+\boldsymbol{\alpha})={[\omega_{1}~{}\cdots~{}\omega_{n}]^{\mathsf{T}}}+\boldsymbol{\beta},$
(14a)

where $\boldsymbol{\beta}$ denotes the correction to the natural frequencies.
Problem (14) always admits a solution because $\boldsymbol{\beta}$ can be
chosen to satisfy the constraint (25a) for any choice of the network
parameters $\boldsymbol{\delta}+\boldsymbol{\alpha}$. Further, the (unique)
solution to the minimization problem (14) can also be computed in closed form:

$\displaystyle\begin{bmatrix}\boldsymbol{\alpha}^{*}\\\
\boldsymbol{\beta}^{*}\end{bmatrix}=\begin{bmatrix}BD({\boldsymbol{x}})&{-}I_{n}\end{bmatrix}^{\dagger}\left({[\omega_{1}~{}\cdots~{}\omega_{n}]^{\mathsf{T}}}-BD({\boldsymbol{x}})\boldsymbol{\delta}\right)$

where $I_{n}$ denotes the $n\times 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
$\ell_{2}$-norm in the cost function (e.g., the $\ell_{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 [60, 59, 14, 45].
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 [82, 53], 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].

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 ${\boldsymbol{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 Figure 7a illustrates the basin of
attraction of $R$, which we numerically estimate to be half of the torus.

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\in\mathbb{R}^{12\times 12}$ of Pearson correlation coefficients (also known
as functional connectivity), and a functional pattern
$R\in\mathbb{R}^{12\times 12}$ (as in equation (2)) from the phases extracted
by solving the _nonconvex phase synchronization_ problem [11]. Strikingly, we
find that $\|F-R\|_{2}\ll 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
$\boldsymbol{\omega}=BD({\boldsymbol{x}})\boldsymbol{\delta}$, where
${\boldsymbol{x}}$ are phase differences obtained from the previous step, and
integrate the Kuramoto model in equation (1) with random initial conditions
close to ${\boldsymbol{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 [83, 47].

### 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 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 structure-preserving power grid model have the same stable
synchronization manifold of equation (1). In this model,
$\omega_{i}=p_{\ell_{i}}/D_{i}$ is the active power load at node $i$, where
$D_{i}$ is the damping coefficient, and
$A_{ij}=|v_{i}||v_{j}|\mathcal{I}({Y}_{ij})/D_{i}$, with $v_{i}$ denoting the
nodal voltage magnitude and $\mathcal{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 where
$D_{i}=1$, which is possibly due to local excitation controllers. Notice that,
when the phase angles $\boldsymbol{\theta}$ are phase-locked and $A_{ij}$ is
fixed, the active power flow is given by $A_{ij}\sin(\theta_{j}-\theta_{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.

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
$\boldsymbol{\theta}$ 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 $\boldsymbol{\delta}$ 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.

We first utilize MATPOWER [89] to solve the power flow problem. Then we use
the active powers $\boldsymbol{p}_{\ell}$ and voltages $v$ at the buses to
obtain the natural frequencies ${\boldsymbol{\omega}}$ and the adjacency
matrix $A$ of the oscillators, respectively, while the voltage phase angles
are used as initial conditions $\boldsymbol{\theta}(0)$ for the Kuramoto model
in equation (1). We integrate the Kuramoto dynamics and let the voltage phases
$\boldsymbol{\theta}(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
${\boldsymbol{x}}_{\mathrm{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 synchronization-based 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.

## Methods

### 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 $B_{k\ell}=-1$ if oscillator $k$ is the source node of the interconnection
$\ell$, $B_{k\ell}=1$ if oscillator $k$ is the sink node of the
interconnection $\ell$, and $B_{k\ell}=0$ otherwise. The matrix form of
equation (1) can be written as

$\displaystyle\dot{\boldsymbol{\theta}}=$
$\displaystyle~{}[\omega_{1}~{}\cdots~{}\omega_{n}]^{\mathsf{T}}-B\begin{bmatrix}\ddots&&\\\
&\sin(x_{ij})&\\\ &&\ddots\end{bmatrix}\boldsymbol{\delta}$ $\displaystyle=$
$\displaystyle~{}[\omega_{1}~{}\cdots~{}\omega_{n}]^{\mathsf{T}}-B\operatorname{diag}(\\{\sin(x_{ij})\\}_{(i,j)\in\mathcal{E}})\boldsymbol{\delta}=[\omega_{1}~{}\cdots~{}\omega_{n}]^{\mathsf{T}}-BD({\boldsymbol{x}})\boldsymbol{\delta},$

where $D({\boldsymbol{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
$\boldsymbol{1}^{\mathsf{T}}B=0$, we have
$\boldsymbol{1}^{\mathsf{T}}\dot{\boldsymbol{\theta}}=\boldsymbol{1}^{\mathsf{T}}k\boldsymbol{1}=\boldsymbol{1}^{\mathsf{T}}[\omega_{1}~{}\cdots~{}\omega_{n}]^{\mathsf{T}}$,
thus showing that, in any phase locked trajectory, the oscillators frequency
$k$ needs to equal the mean natural frequency
$\frac{1}{n}\sum_{i=1}^{n}\omega_{i}$.

### Any feasible functional pattern has $n-1$ degrees of freedom

The 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
$\boldsymbol{x}=M\theta$, where $M\in\mathbb{R}^{|\mathcal{E}|\times n}$ is
the matrix whose $k$-th row, associated to $x_{ij}$, is all zeros except for
$b_{ki}=-1$ and $b_{kj}=1$. Consider the first $n-1$ rows of $M$, associated
to $x_{12},x_{13},\dots x_{1n}$, and notice that they are linearly
independent. Moreover, the row associated to $x_{ij}$, $i>1$, can be obtained
by subtracting the row associated to $x_{1i}$ to the row associated to
$x_{1j}$, 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
$M_{\text{min}}$ (e.g., any $n-1$ rows corresponding to the transpose
incidence matrix of a spanning tree [30]). We let
$\boldsymbol{x}_{\text{min}}=M_{\text{min}}\boldsymbol{\theta}$, where
$\boldsymbol{x}_{\text{min}}$ is a smallest set of phase differences that can
be used to quantify the synchronization angles among all oscillators. Because
$\ker(M_{\text{min}})=\boldsymbol{1}$, we can obtain the phases
$\boldsymbol{\theta}$ from $\boldsymbol{x}_{\text{min}}$ modulo rotation:
$\boldsymbol{\theta}=M_{\text{min}}^{\dagger}\boldsymbol{x}_{\text{min}}-c\boldsymbol{1}$,
where $M_{\text{min}}^{\dagger}$ denotes the Moore-Penrose pseudo-inverse of
$M_{\text{min}}$ and $c$ is some real number. Further, since
$\ker(M_{\text{min}})=\ker(M)$, we can reconstruct all phase differences
$\boldsymbol{x}$ from $\boldsymbol{x}_{\text{min}}$:

$MM_{\text{min}}^{\dagger}\boldsymbol{x}_{\text{min}}=M(\boldsymbol{\theta}+c\boldsymbol{1})=M\boldsymbol{\theta}+0=\boldsymbol{x}.$

The above equation reveals that all the differences $\boldsymbol{x}$ are
encoded in $\boldsymbol{x}_{\text{min}}$. That is, any $x_{ij}$ can be written
as a linear combination of the elements in $\boldsymbol{x}_{\mathrm{min}}$.
For example, if $n=3$ and
$\boldsymbol{x}_{\text{min}}=[x_{12}~{}x_{23}]^{\mathsf{T}}$, then $x_{13}$ is
a linear combination of the differences in $\boldsymbol{x}_{\text{min}}$,
i.e., $\boldsymbol{x}=\left[\begin{smallmatrix}-1&1&0\\\ 0&-1&1\\\
-1&0&1\end{smallmatrix}\right]\left[\begin{smallmatrix}-1&1&0\\\
0&-1&1\end{smallmatrix}\right]^{\dagger}\boldsymbol{x}_{\text{min}}$, in which
$x_{13}=x_{12}+x_{23}$. 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 (5)

Rewrite the pattern assignment problem
$BD({\boldsymbol{x}})\boldsymbol{\delta}={\boldsymbol{\omega}}$ as

$\displaystyle
BD({\boldsymbol{x}})\boldsymbol{\delta}=B_{:,\mathcal{H}}D_{\mathcal{H},\mathcal{H}}({\boldsymbol{x}})\boldsymbol{\delta}_{\mathcal{H}}+B_{:,\tilde{\mathcal{H}}}D_{\tilde{\mathcal{H}},\tilde{\mathcal{H}}}({\boldsymbol{x}})\boldsymbol{\delta}_{\tilde{\mathcal{H}}}=\boldsymbol{\omega},$

where the subscripts $\mathcal{H}$ and $\tilde{\mathcal{H}}$ denote the
entries corresponding to the Hamiltonian path in conditions (ii.a)-(ii.b) and
the remaining ones, respectively. Since
$\mathrm{Im}(\bar{B}_{:,\mathcal{H}})=\mathrm{Im}(B_{:,\mathcal{H}}D_{:,\mathcal{H}}({\boldsymbol{x}}))=\mathrm{span}(\textbf{1})^{\perp}$,
$\mathrm{Im}(B_{:,\tilde{\mathcal{H}}}D_{:,\tilde{\mathcal{H}}}({\boldsymbol{x}}))\subseteq\mathrm{span}(\textbf{1})^{\perp}$,
and $\boldsymbol{\omega}\in\mathrm{span}(\textbf{1})^{\perp}$, for any vector
$\boldsymbol{\delta}_{\tilde{\mathcal{H}}}$, the following set of weights
solves the above equation:

$\displaystyle\boldsymbol{\delta}_{\mathcal{H}}$
$\displaystyle=\left(B_{:,\mathcal{H}}D_{\mathcal{H},\mathcal{H}}({\boldsymbol{x}})\right)^{\dagger}\left(\boldsymbol{\omega}-B_{:,\tilde{\mathcal{H}}}D_{\tilde{\mathcal{H}},\tilde{\mathcal{H}}}({\boldsymbol{x}})\boldsymbol{\delta}_{\tilde{\mathcal{H}}}\right)=(D_{\mathcal{H},\mathcal{H}}({\boldsymbol{x}})B_{:,\mathcal{H}}^{\mathsf{T}}B_{:,\mathcal{H}}D_{\mathcal{H},\mathcal{H}}({\boldsymbol{x}}))^{-1}D_{\mathcal{H},\mathcal{H}}({\boldsymbol{x}})B_{:,\mathcal{H}}^{\mathsf{T}}\left(\boldsymbol{\omega}-B_{:,\tilde{\mathcal{H}}}D_{\tilde{\mathcal{H}},\tilde{\mathcal{H}}}({\boldsymbol{x}})\boldsymbol{\delta}_{\tilde{\mathcal{H}}}\right)$

Because the matrix
$D_{\mathcal{H},\mathcal{H}}({\boldsymbol{x}})B_{:,\mathcal{H}}^{\mathsf{T}}B_{:,\mathcal{H}}D_{\mathcal{H},\mathcal{H}}({\boldsymbol{x}})$
is an M-matrix, its inverse has nonnegative entries. Further, by condition
(ii.b),
$D_{\mathcal{H},\mathcal{H}}({\boldsymbol{x}})B_{:,\mathcal{H}}^{\mathsf{T}}\boldsymbol{\omega}$
is strictly positive. Then, the vector
$(D_{\mathcal{H},\mathcal{H}}({\boldsymbol{x}})B_{:,\mathcal{H}}^{\mathsf{T}}B_{:,\mathcal{H}}D_{\mathcal{H},\mathcal{H}}({\boldsymbol{x}}))^{-1}D_{\mathcal{H},\mathcal{H}}({\boldsymbol{x}})B_{:,\mathcal{H}}^{\mathsf{T}}\boldsymbol{\omega}$
is also strictly positive, and so is the solution vector
$\boldsymbol{\delta}_{\mathcal{H}}$ for any sufficiently small and positive
vector $\delta_{\tilde{\mathcal{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\Delta$, where
$\Delta=\operatorname{diag}(\\{\mathrm{sign}(\cos(x_{ij}))\\}_{(i,j)\in\mathcal{E}})$
is a matrix that changes the sign of the columns of $B$ associated to negative
weights in the cosine-scaled network:

$B\Delta D({\boldsymbol{x}})\boldsymbol{\delta}={\boldsymbol{\omega}}$

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 $\Delta\boldsymbol{\delta}$.

### 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\in\mathbb{C}^{n\times n}$ be a complex
matrix. The $i$-th Gerschgorin disk is $\mathcal{D}_{i}=(M_{ii},r_{i})$,
$i=1,\dots,n$, where the radius is $r_{i}=\sum_{j\neq i}|M_{ij}|$ and the
center is $M_{ii}$.

###### Theorem 2

_( Gerschgorin [29])_. The eigenvalues of the matrix $M$ lie within the union
$\bigcup_{i=1}^{n}D_{i}$ of its Gerschgorin disks.

Whenever all target phase differences in ${\boldsymbol{x}}$ satisfy
$|x_{ij}|\leq\frac{\pi}{2}$, the Gerschgorin disks of the Jacobian in equation
(10) all lie in the closed left half-plane. However, for patterns
${\boldsymbol{x}}$ containing phase differences $|x_{ij}|\geq\frac{\pi}{2}$,
the union of the Gerschgorin disks intersects the right half-plane. Reducing
the magnitude of the entries satisfying $A_{ij}\cos(x_{ij})<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 $A_{ij}\cos(x_{ij})<0$ can be removed, so that all the
Gerschgorin disks lie completely in the left-half plane.

### Data availability

The brain data that support the findings of this study are available from the
Supplementary Information of [60],
https://doi.org/10.1371/journal.pcbi.1004100.s006. The IEEE 39 New England
data parameters and interconnection scheme analyzed in this study for the
structure-preserving power network model can be found in the reference
textbook [64] (see also Supplementary Information for modeling assumptions).
The parameters for the simulations on the network-reduced IEEE 39 New England
test case in the Supplementary Information have been obtained from Ref. [77]
and Ref. [48].

### Code availability

Source code and documentation for the numerical simulations presented here are
freely available in GitHub at:
https://github.com/tommasomenara/functional_control with the identifier
10.5281/zenodo.4546413.

## References

  * [1] Hurricane Ida traps Louisianans, shatters the power grid.  Associated Press, Aug 2021.
  * [2] Juan A Acebrón, Luis L Bonilla, Conrad J Pérez Vicente, Félix Ritort, and Renato Spigler.  The Kuramoto model: A simple paradigm for synchronization phenomena.  Reviews of modern physics, 77(1):137, 2005.
  * [3] S. Achard and E. Bullmore.  Efficiency and cost of economical brain functional networks.  PLOS Computational Biology, 3(2):e17, 2007.
  * [4] A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou.  Synchronization in complex networks.  Physics Reports, 469(3):93–153, 2008.
  * [5] A. Arenas, A. Díaz-Guilera, and C. J. Pérez-Vicente.  Synchronization reveals topological scales in complex networks.  Physical Review Letters, 96:114102, Mar 2006.
  * [6] L. Arola-Fernández, A. Díaz-Guilera, and A. Arenas.  Synchronization invariance under network structural transformations.  Physical Review E, 97:060301, Jun 2018.
  * [7] 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.
  * [8] K. Bansal, J. O. Garcia, S. H. Tompson, T. Verstynen, J. M. Vettel, and S. F. Muldoon.  Cognitive chimera states in human brain networks.  Science Advances, 5(4), 2019.
  * [9] M. Beliaev, D. Zöllner, A. Pacureanu, P. Zaslansky, and I. Zlotnikov.  Dynamics of topological defects and structural synchronization in a forming periodic tissue.  Nature Physics, 2021.
  * [10] I. Belykh, R. Jeter, and V. Belykh.  Foot force models of crowd dynamics on a wobbly bridge.  Science Advances, 3(11):e1701512, 2017.
  * [11] N. Boumal.  Nonconvex phase synchronization.  SIAM Journal on Optimization, 26(4):2355–2377, 2016.
  * [12] J. Buck and E. Buck.  Biology of synchronous flashing of fireflies.  Nature, 211(5049):562–564, 1966.
  * [13] J. W. Busby, K. Baker, M. D. Bazilian, A. Q. Gilbert, E. Grubert, V. Rai, J. D. Rhodes, S. Shidore, C. A. Smith, and M. E. Webber.  Cascading risks: Understanding the 2021 winter blackout in Texas.  Energy Research & Social Science, 77:102106, 2021.
  * [14] J. Cabral, E. Hugues, O. Sporns, and G. Deco.  Role of local network oscillations in resting-state functional connectivity.  NeuroImage, 57(1):130–139, 2011.
  * [15] A. Cavagna, A. Cimarelli, I. Giardina, G. Parisi, R. Santagati, F. Stefanini, and M. Viale.  Scale-free correlations in starling flocks.  Proceedings of the National Academy of Sciences, 107(26):11865–11870, 2010.
  * [16] S. Chen, M. Fazlyab, M. Morari, G. J. Pappas, and V. M. Preciado.  Learning region of attraction for nonlinear systems.  arXiv preprint arXiv:2110.00731, 2021.
  * [17] M. C. Cross and P. C. Hohenberg.  Pattern formation outside of equilibrium.  Reviews of Modern Physics, 65(3):851, 1993.
  * [18] A. Daffertshofer and B. van Wijk.  On the influence of amplitude on the connectivity between phases.  Frontiers in Neuroinformatics, 5:6, 2011.
  * [19] G. Deco, J. Cruzat, J. Cabral, E. Tagliazucchi, H. Laufs, N. K. Logothetis, and M. L. Kringelbach.  Awakening: Predicting external stimulation to force transitions between different brain states.  Proceedings of the National Academy of Sciences, 116(36):18088–18097, 2019.
  * [20] G. Deco, V. Jirsa, A. R. McIntosh, O. Sporns, and R. Kötter.  Key role of coupling, delay, and noise in resting brain fluctuations.  Proceedings of the National Academy of Sciences, 106(25):10302–10307, 2009.
  * [21] F. Dörfler and F. Bullo.  Synchronization and transient stability in power networks and nonuniform Kuramoto oscillators.  SIAM Journal on Control and Optimization, 50(3):1616–1642, 2012\.
  * [22] F. Dörfler and F. Bullo.  Synchronization in complex networks of phase oscillators: A survey.  Automatica, 50(6):1539–1564, 2014.
  * [23] F. Dörfler, 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.
  * [24] R. M. D’Souza and M. Mitzenmacher.  Local cluster aggregation models of explosive percolation.  Physical Review Letters, 104:195702, May 2010.
  * [25] M. Fazlyab, F. Dörfler, and V. M. Preciado.  Optimal network design for synchronization of coupled oscillators.  Automatica, 84:181–189, 2017.
  * [26] Juergen Fell and Nikolai Axmacher.  The role of phase synchronization in memory processes.  Nature Reviews Neuroscience, 12(2):105–118, 2011.
  * [27] D. R. Fields, A. Araque, H. Johansen-Berg, S.-S. Lim, G. Lynch, K.-A. Nave, M. Nedergaard, R. Perez, T. Sejnowski, and H. Wake.  Glial biology in learning and cognition.  The Neuroscientist, 20(5):426–431, 2014.
  * [28] A. Forrow, F. G. Woodhouse, and J. Dunkel.  Functional control of network dynamics using designed Laplacian spectra.  Physical Review X, 8(4):041043, 2018.
  * [29] S. Geršgorin.  Über die abgrenzung der eigenwerte einer matrix.  Bulletin de l’Académie des Sciences de l’URSS. Classe des sciences mathématiques et na, pages 749–754, 1931.
  * [30] C. Godsil and G. F. Royle.  Algebraic Graph Theory.  Graduate Texts in Mathematics. Springer New York, 2001.
  * [31] Martin Golubitsky, Ian Stewart, and Andrei Török.  Patterns of synchrony in coupled cell networks with multiple arrows.  SIAM Journal on Applied Dynamical Systems, 4(1):78–100, 2005.
  * [32] M. Grant, S. Boyd, and Y. Ye.  CVX: Matlab software for disciplined convex programming, 2009.
  * [33] F. Hellmann, P. Schultz, P. Jaros, R. Levchenko, T. Kapitaniak, J. Kurths, and Y. Maistrenko.  Network-induced multistability through lossy coupling and exotic solitary states.  Nature Communications, 11(1):592, 2020.
  * [34] H. Hong and S. H. Strogatz.  Kuramoto model of coupled oscillators with positive and negative coupling parameters: An example of conformist and contrarian oscillators.  Physical Review Letters, 106(5):054102, 2011.
  * [35] F. C. Hoppensteadt and E. M. Izhikevich.  Weakly Connected Neural Networks.  Springer, 1997.
  * [36] R. A. Horn and C. R. Johnson.  Matrix Analysis.  Cambridge University Press, 1985.
  * [37] P. Hövel, A. Viol, P. Loske, L. Merfort, and V. Vuksanović.  Synchronization in functional networks of the human brain.  Journal of Nonlinear Science, pages 1–24, 2018.
  * [38] S. Kaka, M. R. Pufall, W. H. Rippard, T. J. Silva, S. E. Russek, and J. A. Katine.  Mutual phase-locking of microwave spin torque nano-oscillators.  Nature, 437(7057):389–392, 2005.
  * [39] C. Kirst, M. Timme, and D. Battaglia.  Dynamic information routing in complex networks.  Nature Communications, 7(1):11061, 2016.
  * [40] István Z. Kiss, Craig G. Rusin, Hiroshi Kori, and John L. Hudson.  Engineering complex dynamical structures: Sequential patterns and desynchronization.  Science, 316(5833):1886–1889, 2007.
  * [41] G. Korniss, M. A. Novotny, H. Guclu, Z. Toroczkai, and P. A. Rikvold.  Suppressing roughness of virtual times in parallel discrete-event simulations.  Science, 299(5607):677–679, 2003.
  * [42] A. L’Abbate, G. Migliavacca, U. Häger, C. Rehtanz, S. Rüberg, H. Ferreira, G. Fulli, and A. Purvins.  The role of FACTS and HVDC in the future paneuropean transmission system development.  In 9th IET International Conference on AC and DC Power Transmission (ACDC 2010), pages 1–8, 2010.
  * [43] G. S. Medvedev and X. Tang.  Stability of twisted states in the Kuramoto model on Cayley and random graphs.  Journal of Nonlinear Science, 25(6):1169–1208, 2015.
  * [44] D. Mehta, N. S. Daleo, F. Dörfler, 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.
  * [45] T. Menara, G. Baggio, D. S. Bassett, and F. Pasqualetti.  A framework to control functional connectivity in the human brain.  In IEEE Conf. on Decision and Control, pages 4697–4704, Nice, France, December 2019.
  * [46] T. Menara, G. Baggio, D. S. Bassett, and F. Pasqualetti.  Stability conditions for cluster synchronization in networks of heterogeneous Kuramoto oscillators.  IEEE Transactions on Control of Network Systems, 7(1):302 – 314, 2020.
  * [47] T. Menara, G. Lisi, F. Pasqualetti, and A. Cortese.  Brain network dynamics fingerprints are resilient to data heterogeneity.  Journal of Neural Engineering, 18(2):026004, 2021.
  * [48] A. Moeini, I. Kamwa, P. Brunelle, and G. Sybille.  Open data IEEE test systems implemented in SimPowerSystems for education and research in power grid dynamics and control.  In 2015 50th International Universities Power Engineering Conference (UPEC), pages 1–6, 2015.
  * [49] Joon-Young Moon, UnCheol Lee, Stefanie Blain-Moraes, and George A Mashour.  General relationship of global topology, local dynamics, and directionality in large-scale brain networks.  PLOS Computational Biology, 11(4):e1004225, 2015.
  * [50] Adilson E Motter, Seth A Myers, Marian Anghel, and Takashi Nishikawa.  Spontaneous synchrony in power-grid networks.  Nature Physics, 9(3):191–197, 2013.
  * [51] T. Nishikawa and A. E. Motter.  Maximum performance at minimum cost in network synchronization.  Physica D: Nonlinear Phenomena, 224(1-2):77–89, 2006.
  * [52] T. Nishikawa and A. E. Motter.  Network synchronization landscape reveals compensatory structures, quantization, and the positive effect of negative interactions.  Proceedings of the National Academy of Sciences, 107(23):10342–10347, 2010.
  * [53] 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.
  * [54] M. J. Panaggio, M. Ciocanel, L. Lazarus, C. M. Topaz, and B. Xu.  Model reconstruction from temporal data for coupled oscillator networks.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(10):103116, 2019.
  * [55] L. Papadopoulos, J. Z. Kim, J. Kurths, and D. S. Bassett.  Development of structural correlations and synchronization from adaptive rewiring in networks of Kuramoto oscillators.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 27(7):073115, 2017.
  * [56] L. M. Pecora, F. Sorrentino, A. M. Hagerstrom, T. E. Murphy, and R. Roy.  Cluster synchronization and isolated desynchronization in complex networks with symmetries.  Nature Communications, 5(1):1–8, 2014.
  * [57] Louis M. Pecora and Thomas L. Carroll.  Master stability functions for synchronized coupled systems.  Physical Review Letters, 80:2109–2112, Mar 1998.
  * [58] A. Pikovsky, M. Rosenblum, and J. Kurths.  Synchronization: A Universal Concept in Nonlinear Sciences.  Cambridge University Press, 2003.
  * [59] A. Politi and M. Rosenblum.  Equivalence of phase-oscillator and integrate-and-fire models.  Physical Review E, 91(4):042916, 2015.
  * [60] A. Ponce-Alvarez, G. Deco, P. Hagmann, G. L. Romani, D. Mantini, and M. Corbetta.  Resting-state temporal synchronization networks emerge from connectivity topology and heterogeneity.  PLoS Computational Biology, 11(2):1–23, 02 2015.
  * [61] 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.
  * [62] Y. Qin, T. Menara, D. S. Bassett, and F. Pasqualetti.  Phase-amplitude coupling in neuronal oscillator networks.  Physical Review Research, 3(2), June 2021.
  * [63] H. Sakaguchi and Y. Kuramoto.  A soluble active rotater model showing phase transitions via mutual entertainment.  Progress of Theoretical Physics, 76(3):576–581, 1986.
  * [64] P. W. Sauer and M. A. Pai.  Power System Dynamics and Stability.  Prentice Hall, 1998.
  * [65] A. Schnitzler and J. Gross.  Normal and pathological oscillatory communication in the brain.  Nature Reviews Neuroscience, 6(4):285, 2005.
  * [66] K. Sharafutdinov, L. Rydin Gorjão, M. Matthiae, T. Faulwasser, and D. Witthaut.  Rotor-angle versus voltage instability in the third-order model for synchronous generators.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 28(3):033117, 2018.
  * [67] P. S. Skardal and A. Arenas.  Control of coupled oscillator networks with application to microgrid technologies.  Science Advances, 1(7):e1500339, 2015.
  * [68] P. S. Skardal and A. Arenas.  Memory selection and information switching in oscillator networks with higher-order interactions.  Journal of Physics: Complexity, 2(1):015003, dec 2020.
  * [69] P. S. Skardal, D. Taylor, J. Sun, and A. Arenas.  Collective frequency variation in network synchronization and reverse PageRank.  Physical Review E, 93:042314, Apr 2016.
  * [70] F. Sorrentino, M. Di Bernardo, and F. Garofalo.  Synchronizability and synchronization dynamics of weighed and unweighed scale free networks with degree mixing.  International Journal of Bifurcation and Chaos, 17(07):2419–2434, 2007.
  * [71] F. Sorrentino and E. Ott.  Network synchronization of groups.  Physical Review E, 76(5):056114, 2007.
  * [72] I. Stewart, M. Golubitsky, and M. Pivato.  Symmetry groupoids and patterns of synchrony in coupled cell networks.  SIAM Journal on Applied Dynamical Systems, 2(4):609–646, 2003.
  * [73] S. H. Strogatz.  Exploring complex networks.  Nature, 410(6825):268–276, 2001.
  * [74] S. H. Strogatz.  SYNC: The Emerging Science of Spontaneous Order.  Hyperion, 2003.
  * [75] S. H. Strogatz and I. Stewart.  Coupled oscillators and biological synchronization.  Scientific American, 269(6):102–109, 1993.
  * [76] J. Sun, E. M. Bollt, and T. Nishikawa.  Master stability functions for coupled nearly identical dynamical systems.  EPL (Europhysics Letters), 85(6):60011, 2009.
  * [77] Y. Susuki, I. Mezić, and T. Hikihara.  Coherent swing instability of power grids.  Journal of nonlinear science, 21(3):403–439, 2011.
  * [78] A. Townsend, M. Stillman, and S. H. Strogatz.  Dense networks that do not synchronize and sparse ones that do.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(8):083142, 2020.
  * [79] K. E. Trenberth.  Spatial and temporal variations of the Southern Oscillation.  Quarterly Journal of the Royal Meteorological Society, 102(433):639–653, 1976.
  * [80] M. P. Van Den Heuvel and H. E. Hulshoff Pol.  Exploring the brain network: a review on resting-state fmri functional connectivity.  European Neuropsychopharmacology, 20(8):519–534, 2010.
  * [81] M. van Wingerden, M. Vinck, J. Lankelma, and C. M. A. Pennartz.  Theta-band phase locking of orbitofrontal neurons during reward expectancy.  Journal of Neuroscience, 30(20):7078–7087, 2010.
  * [82] F. Varela, J. P. Lachaux, E. Rodriguez, and J. Martinerie.  The brainweb: Phase synchronization and large-scale integration.  Nature Reviews Neuroscience, 2(4):229–239, 2001.
  * [83] F. Váša, M. Shanahan, P. J. Hellyer, G. Scott, J. Cabral, and R. Leech.  Effects of lesions on synchrony and metastability in cortical networks.  Neuroimage, 118:456–467, 2015.
  * [84] J.-W. Wang and L.-L. Rong.  Cascade-based attack vulnerability on the US power grid.  Safety Science, 47(10):1332 – 1336, 2009.
  * [85] D. Witthaut and M. Timme.  Braess’s paradox in oscillator networks, desynchronization and power outage.  New Journal of Physics, 14(8):083036, aug 2012.
  * [86] D. Zelazo and M. Bürger.  On the robustness of uncertain consensus networks.  IEEE Transactions on Control of Network Systems, 4(2):170–178, 2017\.
  * [87] F. Zhang and N. E. Leonard.  Coordinated patterns of unit speed particles on a closed curve.  Systems & Control Letters, 56(6):397–407, 2007.
  * [88] L. Zhu and D. J. Hill.  Synchronization of Kuramoto oscillators: A regional stability framework.  IEEE Transactions on Automatic Control, 65(12):5070–5082, 2020\.
  * [89] R. D. Zimmerman, C. E. Murillo-Sánchez, 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.

## Supplementary Information

##  1 Supplementary Text

###  1.1 Comparison between our correlation metric and the Pearson correlation
coefficient

In this work, we provide methods and principles to guarantee the emergence of
phase-locked trajectories that describe functional patterns. The latter are
defined utilizing $\rho_{ij}=<\cos(\theta_{j}-\theta_{i})>_{t}$, instead of
the classical Pearson correlation coefficient. The same metric is also
utilized in [19] to quantify the similarity between BOLD signals in functional
MRI recordings, and in [5], to quantify the synchrony in hierarchical
networks.

To compare the two metrics, recall the definition of the classical (sample)
Pearson correlation, which is a measure of linear correlation between two
series of data, for vectors of samples $y\in\mathbb{R}^{T}$ and
$z\in\mathbb{R}^{T}$:

$r_{ij}=\frac{\sum_{i=1}^{N}(y_{i}-y_{\mathrm{mean}})(z_{i}-z_{\mathrm{mean}})}{\sqrt{\sum_{i=1}^{N}(y_{i}-y_{\mathrm{mean}})^{2}}\sqrt{\sum_{i=1}^{N}(z_{i}-z_{\mathrm{mean}})^{2}}},$

where $y_{\mathrm{mean}}$ and $z_{\mathrm{mean}}$ denote the sample means of
vectors $y$ and $z$, respectively. Notice that the length $T$ of the vectors
$y$ and $z$ depends on the sampling time and on the window length. We argue
that our metric (equation (2) in the main text) is simply more convenient when
dealing with periodic phase signals. While we could use the Pearson
correlation coefficient to define functional pattern instead of our cosine-
based metric, the former does not perform well on periodic signals that evolve
on the unit circle, and is heavily dependent on the time window employed to
collect the samples. Thus, specific adjustments are needed to define
correlation patterns through Pearson correlation coefficient for the class of
time series (phase trajectories) studied in this work.

As an example, consider two identical sinusoidal signals (i.e., phase-locked)
with natural frequency $\omega=2\pi$ but shifted initial conditions:
$\theta_{1}(0)=0$, $\theta_{2}(0)=\varphi$, with $\varphi\in(0,~{}\pi]$. Fig.
10 illustrates the differences between the Pearson correlation coefficient
$r_{12}$ and our correlation metric $\rho_{12}$ computed over a time window of
varying length and for different values of the initial phase shift $\varphi$.
In all panels, the values of $r_{12}$ vary at each point, emphasizing the
dependence of the Pearson correlation coefficient from the length of the time
window. Conversely, in all panels, the value of $\rho_{12}$ remains unaltered
by the choice of time window length.

In conclusion, we choose to define functional patterns through
$\rho_{ij}=<\cos(\theta_{j}-\theta_{i})>_{t}$ because it is a convenient and
suitable metric for this class of phase trajectories that naturally emerge in
oscillator systems.

###  1.2 Stability results for functional patterns with angle differences
larger than $\frac{\pi}{2}$ in the case of line and cycle with positive
weights

Consider a line network of $n$ (ordered) oscillators with positive-only
weights that possesses an equilibrium for the phase difference dynamics
satisfying $|x_{i,i+1}|>\frac{\pi}{2}$ for some $i\in\\{1,\dots,n-1\\}$. It is
straightforward to deduce from the results on structural balance [86] that the
considered equilibrium is unstable. This implies that the functional 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}|>\frac{\pi}{2}$ and allowing
only positive weights is the cycle (i.e., a line network where the first and
last oscillators are connected). Consider a cycle network of $n>4$ oscillators
with positive weights, and denote with ${\boldsymbol{x}}_{\mathrm{cycle}}$ a
vector of the $n$ phase differences between connected oscillators. We find
that, after a suitable relabeling of the oscillators for which $x_{12}$
satisfies $|x_{12}|>\frac{\pi}{2}$:

###### Theorem 2

(Stability of phase differences equilibria with $|x_{ij}|>\frac{\pi}{2}$ in
cycle networks) The equilibrium
${\boldsymbol{x}}_{\mathrm{cycle}}=[x_{12}~{}x_{23}~{}\dots~{}x_{n-1,n}]^{\mathsf{T}}=[\gamma~{}\varphi_{1}~{}\dots~{}\varphi_{n-1}]^{\mathsf{T}}$
with $|\gamma|>\frac{\pi}{2}$, $|\varphi_{i}|<\frac{\pi}{2}$, and
$\varphi_{i}=\varphi_{n-i}$ for all $i=1,\dots,n-1$, is stable if and only if

  * _(i)_

$A_{i,i+1}=A_{12}\frac{\sin(\gamma)}{\sin(\varphi_{i-1})}$
for $i=2,\dots,n$, with $n-i+3\triangleq 1$ if $i=2$, and ${n+1}\triangleq 1$;

  * _(ii)_

$|\cot{\gamma}|\leq\left(\tan(\varphi_{1})+\dots+\tan(\varphi_{n-1})\right)^{-1}$.

Moreover, if $\varphi_{1}=\dots=\varphi_{n-1}$ and $n\to\infty$, the largest
possible value for $\gamma$ such that $\boldsymbol{x}_{\mathrm{cycle}}$ is
stable tends to the value $\gamma\approx 1.789776$, which is the solution to
$\gamma-\tan(\gamma)=2\pi$.

Proof of Theorem 2: To assess the stability of the equilibrium
$\bar{\boldsymbol{x}}_{\mathrm{cycle}}=[x_{12}~{}x_{23}~{}\dots~{}x_{n-1,n}]^{\mathsf{T}}=[\gamma~{}\varphi_{1}~{}\dots~{}\varphi_{n-1}]^{\mathsf{T}}$
with $|\gamma|>\frac{\pi}{2}$, $|\varphi_{i}|<\frac{\pi}{2}$, and
$\varphi_{i}=\varphi_{n-i}$ for all $i=1,\dots,n-1$, we analyze the spectrum
of the Jacobian
$J({\boldsymbol{x}}_{\mathrm{cycle}})=-\mathcal{L}({\boldsymbol{x}}_{\mathrm{cycle}})$.
From Ref. [86, Corollary IV.7], we have that a necessary and sufficient
condition for the Laplacian matrix
$\mathcal{L}({\boldsymbol{x}}_{\mathrm{cycle}})$ of the cosine-scaled network
to be positive semidefinite is

$|{A_{12}}\cos(\gamma)|\leq\mathcal{R}_{12}^{-1},$ (15)

with $\mathcal{R}_{12}$ being the effective resistance of the graph in which
the edge $(1,2)$ has been removed. That is,

$\mathcal{R}_{12}=\frac{1}{{A_{23}}\cos(\varphi_{1})}+\dots+\frac{1}{{A_{n-1,n}}\cos(\varphi_{n-1})}.$
(16)

Since the adjacency matrix satisfies $A=A^{\mathsf{T}}$ and
${\boldsymbol{x}}_{\mathrm{cycle}}$ is an equilibrium for the difference
dynamics of the cycle network, the network weights must be identical pairwise:

$A_{i,i+1}=A_{n-i+2,n-i+3}=A_{12}\frac{\sin(\gamma)}{\sin(\varphi_{i-1})},$
(17)

for $i=2,\dots,n$ with the convention $n-i+3\triangleq 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

$\mathcal{R}_{12}=\frac{\tan(\varphi_{1})}{{A_{12}}\sin(\gamma)}+\dots+\frac{\tan(\varphi_{n-1})}{{A_{12}}\sin(\gamma)},$
(18)

which makes the condition in Eq. (15) become, after algebraic calculations,
condition (ii) of Theorem 2:

$|\cot{\gamma}|\leq\left(\tan(\varphi_{1})+\dots+\tan(\varphi_{n-1})\right)^{-1}.$
(19)

Thus, given that the Laplacian
$\mathcal{L}({\boldsymbol{x}}_{\mathrm{cycle}})$ is positive definite, the
Jacobian $J({\boldsymbol{x}}_{\mathrm{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 $\varphi_{1}=\dots=\varphi_{n-1}$, we have that

$\varphi_{i}=\frac{2\pi-\gamma}{n-1}.$

Hence, the right-hand side of equation (19) becomes
$\frac{\cot(\varphi)}{n-1}$, and
$\lim_{n\to\infty}\frac{\cot(\varphi)}{n-1}=\frac{1}{2\pi-\gamma}$. Since
$|\gamma|>\frac{\pi}{2}$, plugging the limit value for $\frac{1}{2\pi-\gamma}$
into equation (19) and solving for the equality yields
$\gamma-\tan(\gamma)=2\pi$, whose unique solution is $\gamma\approx 1.789776$.
This concludes the proof. $acksquare$

###  1.3 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}\geq 0$, there is no unique pattern. In fact, the largest basin of
attraction is achieved by $\mathbb{S}^{1}$_-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
${\boldsymbol{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: $\omega_{1}=\omega_{2}$. The only stable equilibrium
for the phase differences dynamics is $x_{12}=\theta_{2}-\theta_{1}=0$. Yet,
there also exists an unstable equilibrium $x_{12}=\pi$. 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 $\mathbb{S}^{1}$_-synchronizing_ clusters with possibly different natural
frequencies. We consider cluster synchronization in networks with a partition
of the oscillators $\mathcal{C}=\\{\mathcal{C}_{1},\dots,\mathcal{C}_{m}\\}$,
with $\mathcal{C}_{k}\subseteq\mathcal{O}$ being a subset of the network
oscillators constituting an $\mathbb{S}^{1}$_-synchronizing_ graph,
$k\in\\{1,\dots,m\\}$, where $\bigcup_{k=1}^{m}\mathcal{C}_{k}=\mathcal{O}$
and $\mathcal{C}_{k}\cap\mathcal{C}_{\ell}=\emptyset$ if $k\neq\ell$. Further,
$\omega_{i}=\omega_{j}$ for every $i,j\in\mathcal{C}_{k}$,
$k\in\\{1,\dots,m\\}$.111This is a necessary and sufficient condition for
cluster synchronization [46]. Whenever cluster synchronization emerges, the
diagonal blocks (of sizes $|\mathcal{C}_{k}|\times|\mathcal{C}_{k}|$,
$k=1,\dots,m$) of the functional pattern associated with partition
$\mathcal{C}$ satisfy $\rho_{ij}=1$ (see Supplementary Fig. 15).222Without
loss of generality, we assume that the network oscillators are labeled such
that consecutive labels belong to the same cluster. Thus, for a given number
of clusters $m\geq 1$, each of the ${n\choose 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
$\mathbb{S}^{1}$_-synchronizing_ graphs, we find that, in networks where at
least one cluster satisfies $|\mathcal{C}_{k}|\geq 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 $\mathbb{S}^{1}$_-synchronizing_ graphs (i.e.,
$x_{ij}=\pi$ 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}=\pi$.

As an example, consider a $4$-oscillator network partitioned as
$\mathcal{C}=\\{\mathcal{C}_{1},\mathcal{C}_{2}\\}$, where
$\mathcal{C}_{1}=\\{1,2\\}$ and $\mathcal{C}_{2}=\\{3,4\\}$. The network
parameters read

$A=\begin{bmatrix}0&5&2&0\\\ 5&0&0&2\\\ 2&0&0&6\\\ 0&2&6&0\end{bmatrix}\
\text{ and }\ \boldsymbol{\omega}=\begin{bmatrix}1\\\ 1\\\ 2\\\ 2\\\
\end{bmatrix}.$

Note that the two clusters are $\mathbb{S}^{1}$_-synchronizing_ graphs. It can
be shown that there exists an unstable equilibrium at
${\boldsymbol{x}}_{\mathrm{desired}}=[x_{12}~{}x_{23}~{}x_{34}]^{\mathsf{T}}=[\pi~{}0.25268~{}\pi]^{\mathsf{T}}$.
Thus, whenever the pattern $R$ associated to the partition $\mathcal{C}$ is
stable, it is at most almost-globally stable.

###  1.4 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
${\boldsymbol{x}}$. These equilibria correspond to phase trajectories $\theta$
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 $\theta_{1}$ as a reference, we choose two points for
the phase differences $x_{1i}=\theta_{i}-\theta_{1}$ to be set as equilibria:
${\boldsymbol{x}}^{(1)}_{\mathrm{desired}}=\left[\frac{\pi}{6}~{}\frac{\pi}{6}~{}\frac{\pi}{4}~{}\frac{\pi}{4}~{}\frac{\pi}{6}~{}\frac{\pi}{4}\right]^{\mathsf{T}}$
and
${\boldsymbol{x}}^{(2)}_{\mathrm{desired}}=\left[\frac{\pi}{8}~{}\frac{\pi}{3}~{}\frac{\pi}{4}~{}\frac{\pi}{4}~{}\frac{\pi}{6}~{}\frac{\pi}{4}\right]^{\mathsf{T}}$.
The initial network parameters (adjacency matrix and zero-mean natural
frequencies) read as:

$A=\begin{bmatrix}0&2&2&2&2&2&2\\\ 2&0&2&2&2&2&2\\\ 2&2&0&2&2&2&2\\\
2&2&2&0&2&2&2\\\ 2&2&2&2&0&2&2\\\ 2&2&2&2&2&0&2\\\ 2&2&2&2&2&2&0\end{bmatrix}\
\text{and}\ {\boldsymbol{\omega}}=\begin{bmatrix}0.3160\\\ -0.1266\\\
0.1437\\\ 0.2771\\\ -0.3593\\\ -0.3363\\\ 0.0854\end{bmatrix},$

respectively. To impose the desired equilibria, we numerically solve the
following convex problem through standard cvx routines [32]:

$\displaystyle\min_{\boldsymbol{\alpha}}$ $\displaystyle\ \ \
\|\boldsymbol{\delta}+\boldsymbol{\alpha}\|_{1}$ subject to $\displaystyle\ \
\begin{bmatrix}BD\left({\boldsymbol{x}}^{(1)}\right)\\\\[5.0pt]
BD\left({\boldsymbol{x}}^{(2)}\right)\end{bmatrix}(\boldsymbol{\delta}+\boldsymbol{\alpha})=\begin{bmatrix}{\boldsymbol{\omega}}\\\
{\boldsymbol{\omega}}\end{bmatrix}.$

The solution $\boldsymbol{\alpha}^{*}$ yields the following corrected
adjacency matrix:

$A_{\mathrm{c}}=\begin{bmatrix}0&4.2841&1.3731&-1.6720&-2.5576&2.4333&-1.9382\\\
4.2841&0&-3.0379&2.9361&2.9041&2&2.9252\\\
1.3731&-3.0379&0&0.7026&0.6949&0.5221&0.7\\\ -1.6720&2.9361&0.7026&0&2&2&2\\\
-2.5576&2.9041&0.6949&2&0&2&2\\\ 2.4333&2&0.5221&2&2&0&2\\\
-1.9382&2.9252&0.7&2&2&2&0\end{bmatrix}.$

An investigation of the Jacobian spectrum (see main text for results on
stability) of the phase differences computed at the two equilibria
${\boldsymbol{x}}^{(1)}_{\mathrm{desired}}$ and
${\boldsymbol{x}}^{(2)}_{\mathrm{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
${\boldsymbol{x}}^{(1)}_{\mathrm{desired}}$ at time $t=0$, and converge to
${\boldsymbol{x}}^{(2)}_{\mathrm{desired}}$ after a perturbation is applied at
$t=50$ to force them out from the equilibrium
${\boldsymbol{x}}^{(1)}_{\mathrm{desired}}$.

###  1.5 A heuristic method to promote stability of functional patterns in
positive networks

To promote stability of functional patterns with phase differences
$|x_{ij}|>\frac{\pi}{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
$\mathcal{L}({\boldsymbol{x}}_{\mathrm{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, $|\mathcal{E}|=9$
interconnections, and desired minimum vector of phase differences

${\boldsymbol{x}}_{\text{desired}}^{(1)}=\left[\frac{21\pi}{32}~{}\frac{\pi}{6}~{}\frac{\pi}{6}~{}\frac{\pi}{8}~{}\frac{\pi}{8}~{}\frac{\pi}{3}\right]^{\mathsf{T}},$

where $x_{ij}=\theta_{j}-\theta_{1}$, $j=2,\dots,7$. Notice that the first
difference $x_{12}>\pi/2$, hence $\cos(x_{12})<0$. Consider the oscillator
network with structural parameters that read as:

$A=\begin{bmatrix}0&0.1706&0&0&0.5796&0&0\\\ 0.1706&0&1.3434&0&0&0&0\\\
0&1.3434&0&2&2&2.2140&0\\\ 0&0&2&0&0&1.2392&0.3432\\\ 0.5796&0&2&0&0&2&0\\\
0&0&2.2140&1.2392&2&0&0\\\ 0&0&0&0.3432&0&0&0\\\ \end{bmatrix}\ \text{and}\
\bar{\boldsymbol{\omega}}=\begin{bmatrix}-0.3424\\\ 0.4683\\\ 0.2023\\\
-0.0099\\\ -0.0393\\\ -0.4507\\\ 0.1716\end{bmatrix}.$ (20)

Such a network admits the following phase-locked equilibrium

${\boldsymbol{x}}_{\text{desired}}^{(0)}=\left[\frac{\pi}{4}~{}\frac{\pi}{6}~{}\frac{\pi}{6}~{}\frac{\pi}{8}~{}\frac{\pi}{8}~{}\frac{\pi}{3}\right]^{\mathsf{T}},$

which differs from ${\boldsymbol{x}}_{\text{desired}}^{(1)}$ only in $x_{12}$,
and has all $x_{ij}$ satisfying $|x_{ij}|<\pi/2$ (thus generating a stable
functional pattern). Supplementary Fig. 13a-b illustrate the stability of
${\boldsymbol{x}}_{\text{desired}}^{(0)}$ 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
${\boldsymbol{x}}_{\text{desired}}^{(1)}$ implies that $\mathcal{N}=\\{1\\}$,
so that $\boldsymbol{\alpha}_{\mathcal{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])

$\displaystyle\min_{\boldsymbol{\alpha}}$ $\displaystyle\ \
\left\|\boldsymbol{\delta}_{\mathcal{N}}+\boldsymbol{\alpha}_{\mathcal{N}}\right\|_{1}$
(21) $\displaystyle\mathrm{subject~{}to}$ $\displaystyle\ \
BD({\boldsymbol{x}})(\boldsymbol{\delta}+\boldsymbol{\alpha})={\boldsymbol{\omega}},$
$\displaystyle\mathrm{and}$ $\displaystyle\ \
(\boldsymbol{\delta}+\boldsymbol{\alpha})\geq 0.$

The optimal $\boldsymbol{\alpha}^{*}$ reads

$\boldsymbol{\alpha}^{*}=[-0.1706~{}0.3152~{}-0.8748~{}59~{}0.9242~{}0~{}0~{}0~{}59]^{\mathsf{T}},$

and the adjusted network adjacency matrix becomes:

$\tilde{A}=\begin{bmatrix}0&{\color[rgb]{1,0,0}\boldsymbol{0}}&0&0&0.8948&0&0\\\
{\color[rgb]{1,0,0}\boldsymbol{0}}&0&{0.4686}&0&0&0&0\\\
0&{0.4686}&0&61&2.9242&2.2140&0\\\ 0&0&61&0&0&1.2392&0.3432\\\
0.8948&0&2.9242&0&0&61&0\\\ 0&0&2.2140&1.2392&61&0&0\\\ 0&0&0&0.3432&0&0&0\\\
\end{bmatrix},$ (22)

where the optimal network correction $\boldsymbol{\alpha}^{*}$ has
disconnected the entry $A_{12}$ (highlighted in red). The matrix $\tilde{A}$
guarantees stability of the functional pattern associated with
${\boldsymbol{x}}_{\text{desired}}^{(1)}$, as we illustrate in Supplementary
Fig. 13c-d.

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
$\mathcal{P}$ (resp., $\mathcal{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:

$\displaystyle\min_{\boldsymbol{\alpha}}$ $\displaystyle\ \
c_{1}\left\|\boldsymbol{\alpha}_{\mathcal{P}}\right\|_{\star}+c_{2}\left\|\boldsymbol{\delta}_{\mathcal{N}}+\boldsymbol{\alpha}_{\mathcal{N}}\right\|_{\star}$
(23) $\displaystyle\mathrm{subject~{}to}$ $\displaystyle\ \
BD({\boldsymbol{x}})(\boldsymbol{\delta}+\boldsymbol{\alpha})={\boldsymbol{\omega}},$
$\displaystyle\ \ (\boldsymbol{\delta}+\boldsymbol{\alpha})\geq 0,$

where $\|\cdot\|_{\star}$ is a desired vector norm, $c_{1},c_{2}>0$ are
arbitrary penalty coefficients, $\boldsymbol{\alpha}_{\mathcal{P}}$ denotes
the entries of the tuning vector $\boldsymbol{\alpha}$ that are associated to
positive weights in the cosine-scaled network, and
$\boldsymbol{\alpha}_{\mathcal{N}}$ denotes the entries of the tuning vector
$\boldsymbol{\alpha}$ that are associated to negative weights
$\boldsymbol{\delta}_{\mathcal{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 $\ell_{1}$-norm. The optimal
$\boldsymbol{\alpha}^{*}$ reads

$\boldsymbol{\alpha}^{*}=[-0.1706~{}0.3152~{}-0.8748~{}0~{}0.9242~{}0~{}0~{}0~{}0]^{\mathsf{T}},$

and the adjusted network adjacency matrix becomes:

$\tilde{A}=\begin{bmatrix}0&{\color[rgb]{1,0,0}\boldsymbol{0}}&0&0&0.8948&0&0\\\
{\color[rgb]{1,0,0}\boldsymbol{0}}&0&0.4686&0&0&0&0\\\
0&0.4686&0&2&2.9242&2.2140&0\\\ 0&0&2&0&0&1.2392&0.3432\\\
0.8948&0&2.9242&0&0&2&0\\\ 0&0&2.2140&1.2392&2&0&0\\\ 0&0&0&0.3432&0&0&0\\\
\end{bmatrix},$

where, as in the procedure above, the optimal network correction $\alpha$ 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 $\boldsymbol{\alpha}^{*}$ is gradually applied. The main differences
with respect to the minimization in (21) is that the norm of
$\boldsymbol{\alpha}^{*}$ is smaller, and the eigenvalues of $\tilde{A}$ are
closer to the original ones in $A$.

###  1.6 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 Kuramoto dynamics in equation (3) of the main text
does not hold anymore and requires a rewriting. In what follows, we use the
subscript “$\mathrm{d}$” to indicate notation associated to _directed_ graphs.
Specifically, $\mathcal{E}_{\mathrm{d}}$ denote the oriented edge set, so that
$D_{\mathrm{d}}\in\mathbb{R}^{|\mathcal{E}_{\mathrm{d}}|\times\mathcal{|}\mathcal{E}_{\mathrm{d}}|}$
and
$\boldsymbol{\delta}_{\mathrm{d}}\in\mathbb{R}^{|\mathcal{E}_{\mathrm{d}}|}$
denote the diagonal matrix of all $\sin(x_{ij})$ with
$(j,i)\in\mathcal{E}_{\mathrm{d}}$, and the vector of all the network weights
$A(i,j)\neq 0$, respectively. Further, let us define
$B_{\mathrm{source}}\in\mathbb{R}^{n\times|\mathcal{E}_{\mathrm{d}}|}$ the
modified incidence matrix whose columns have nonzero entries only at the
edges’ sources. That is, $B_{\mathrm{source},{k\ell}}=-1$ if $k$ is the source
of the interconnection $\ell$, and $B_{\mathrm{source},{k\ell}}=0$ otherwise.
These definitions allow us to define the matrix form for a directed network of
Kuramoto oscillators:

$\dot{\boldsymbol{\theta}}=[\omega_{1}~{}\cdots~{}\omega_{n}]^{\mathsf{T}}-B_{\mathrm{source}}D_{\mathrm{d}}(\boldsymbol{x})\boldsymbol{\delta}_{\mathrm{d}}.$
(24)

The main change in the behavior of directed networks with respect to
undirected ones is that the frequencies $\dot{\boldsymbol{\theta}}$ of phase-
locked trajectories do not typically converge to the average natural frequency
(i.e., $\dot{\boldsymbol{\theta}}\neq\omega_{\mathrm{mean}}\boldsymbol{1}$).
For phase-locked trajectories we have that
$\dot{\boldsymbol{\theta}}=\omega_{\mathrm{sync}}\boldsymbol{1}$, where the
constant $\omega_{\mathrm{sync}}\in\mathbb{R}$ is not known _a priori_ , and
can only be estimated in the almost-fully synchronized regime
$|\theta_{i}(t)-\theta_{j}(t)|\ll 1$ for all $i,j\in\mathcal{O}$ [69]. Yet,
not knowing $\omega_{\mathrm{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 $\bar{\boldsymbol{x}}_{\mathrm{desired}}$
and assign a desired phase-locked frequency $\omega_{\mathrm{sync}}$. To do
so, we utilize equation (24) above with
$\dot{\boldsymbol{\theta}}=\omega_{\mathrm{sync}}\boldsymbol{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

$\displaystyle\min_{\boldsymbol{\alpha}}$
$\displaystyle\left\|\boldsymbol{\alpha}\right\|_{1}$ (25)
$\displaystyle\mathrm{subject~{}to}$
$\displaystyle\boldsymbol{\omega}-B_{\mathrm{source}}D_{\mathrm{d}}(\boldsymbol{x})(\boldsymbol{\delta}_{\mathrm{d}}+\boldsymbol{\alpha})=\omega_{\mathrm{sync}}\boldsymbol{1},$
(25a)

for a network of $n=7$ oscillators with adjacency matrix and natural
frequencies as follows:

$A_{\mathrm{d}}=\begin{bmatrix}0&{\color[rgb]{1,0,0}\boldsymbol{0}}&1&1&1&1&1\\\
1&0&1&1&1&1&1\\\ 1&1&0&1&1&1&1\\\ 1&1&1&0&1&1&1\\\ 1&1&1&1&0&1&1\\\
1&1&1&1&1&0&1\\\ 1&1&1&1&1&1&0\end{bmatrix}\ \text{and}\
{\boldsymbol{\omega}}=0.1\cdot\begin{bmatrix}1\\\ 2\\\ 3\\\ 4\\\ 5\\\ 6\\\
7\end{bmatrix},$

where the entry highlighted in red is the one causing the asymmetry in the
network coupling. The target phase differences $x_{1i}$ are
${\boldsymbol{x}}_{\mathrm{desired}}=\left[\frac{\pi}{3}~{}\frac{\pi}{4}~{}\frac{\pi}{6}~{}\frac{\pi}{8}~{}\frac{\pi}{8}~{}\frac{\pi}{6}\right]^{\mathsf{T}}$,
and the desired phase locking frequency is $\omega_{\mathrm{sync}}=1$ rad/sec.
The solution $\boldsymbol{\alpha}^{*}$ to problem (25) above yields

$A_{\mathrm{d}}=\begin{bmatrix}0&0&-1.2238&1&1&1&1\\\ -3.7832&0&1&1&1&1&1\\\
-2.4384&1&0&1&1&1&1\\\ 0.3978&1.6022&1&0&1&1&1\\\ 1&0.3925&1&1&0&1&1\\\
1&0.2282&1&1&1&0&1\\\ 0.6978&1.3022&1&1&1&1&0\end{bmatrix}.$

Supplementary Fig. 12 illustrates that the phase trajectories associated with
such a solution achieve the target functional pattern.

We conclude this discussion by remarking that, besides our optimization
problems, the sufficient conditions (i.a), (i.b), and (i.c) in the main text
for the existence of positive coupling strengths that realize a target pattern
can also be adapted to the case of directed networks. To show this, we let
${\boldsymbol{\omega}}_{\mathrm{d}}=[\omega_{1}-\omega_{\text{sync}}~{}\cdots~{}\omega_{n}-\omega_{\text{sync}}]^{\mathsf{T}}$
and
$\bar{B}_{\mathrm{source}}=B_{\mathrm{source}}~{}\mathrm{sign}(D_{\mathrm{d}}(\boldsymbol{x}))$.
Then, a sufficient condition for the existence of positive network weights
that achieve a desired functional pattern is the following one.

_There exists ${\boldsymbol{\delta}}\geq 0$ such that
$B_{\mathrm{source}}D_{\mathrm{d}}({\boldsymbol{x}}){\boldsymbol{\delta}_{\mathrm{d}}}={\boldsymbol{\omega}}_{\mathrm{d}}$
if there exists a set $\mathcal{S}$ satisfying:_

  1. (iii.a)

_$D_{\mathrm{d}_{ii}}(\boldsymbol{x})D_{\mathrm{d}_{jj}}(\boldsymbol{x})B_{\mathrm{source}_{:,i}}^{\mathsf{T}}B_{\mathrm{source}_{:,j}}\leq
0$ for all $i,j\in\mathcal{S}$ with $i\neq j$;_

  2. (iii.b)

_${\boldsymbol{\omega}}_{\mathrm{d}}^{\mathsf{T}}B_{\mathrm{source}_{:,i}}D_{\mathrm{d}_{ii}}(\boldsymbol{x})
>0$ for all $i\in\mathcal{S}$;_

  3. (iii.c)

_${\boldsymbol{\omega}}_{\mathrm{d}}\in\mathrm{Im}(B_{\mathrm{source}_{:,\mathcal{S}}})$._

Note that, since $\omega_{\mathrm{sync}}$ is not known _a priori_ in most
cases, these conditions are hard to check.

###  1.7 Coupled Kuramoto oscillators to approximate fMRI data

The interaction between static large-scale structural architecture of the
human brain and local oscillations of neural communities is a key factor in
the functional connectivity patterns that are empirically observed through
functional magnetic resonance imaging (fMRI) when the brain is in a resting-
state condition [80]. In the last two decades, extensive literature has
resorted to Kuramoto phase oscillators to model fMRI data [20, 60, 59, 14, 83,
37, 45]. Many works, such as Ref. [60] and Ref. [14], focus on the analysis of
the oscillatory behaviors of neural populations that lead the emergence of
functionally connected networks by modeling fMRI data as the output of
networks of Kuramoto oscillators. The main working assumption is that at each
node of the structural brain network there exists a community of excitatory
and inhibitory neurons whose dynamical state is in a regime of self-sustained
oscillations. From a modeling standpoint, this assumption is equivalent to
employing a network of weakly coupled Wilson-Cowan oscillators [35, 18], or to
a supercritical Andronov-Hopf bifurcation, such as the Stuart-Landau model in
oscillatory regime [49]. In this setting, the neurons’ firing rates describe a
closed periodic trajectory in phase space; that is, the firing rates delineate
a limit cycle. Thus, the dynamics can be approximated by a single variable,
which is the angle (or _phase_) on this cycle. This regime is then modeled by
a network of coupled heterogeneous Kuramoto oscillators that are connected to
each other according to the architecture of the human brain.

###  1.8 Procedure to extract phase-locked trajectories from fMRI data

Following the procedure in Ref. [60], we apply a narrow-band filter in the low
frequency range $[0.04~{}0.07]$ Hz to the time series for each brain region.
Next, to obtain a measure of the functional synchrony between the brain
regions, we generate the $n\times n$ functional connectivity matrix ${F}$,
whose entry ${F}_{ij}$ indicates the pairwise Pearson correlation coefficient
between filtered time series of recorded neural activity. To map these
functional correlations to the phase domain, we extracted the phase time
series $\tilde{\boldsymbol{\theta}}(t)$ by applying a Hilbert Transform to the
filtered signals. From $\tilde{\boldsymbol{\theta}}(t)$, one can identify time
windows over which frequency synchronization emerges with the aid of the
phase-locking value matrix $P=[P_{ij}]$, where

$P_{ij}(t_{0},t_{f})=\frac{1}{t_{f}-t_{0}}\sum_{t=t_{0}}^{t_{f}}\left|e^{\mathsf{i}[\theta_{j}(t)-\theta_{i}(t)]}\right|.$

Clearly, if $P_{ij}(t_{0},t_{f})\approx 1$ for all $i,j$, then the time window
$[t_{0}~{}t_{f}]$ comprises phase-locked trajectories.

Since the phase time series $\tilde{\boldsymbol{\theta}}(t)$ are derived from
inherently noisy measurements, we compute the best estimate of the phases
${\boldsymbol{\theta}}^{*}$ (modulo rotation) that are compatible with the
noisy measurement in $\tilde{\boldsymbol{\theta}}(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 $\tilde{\boldsymbol{\theta}}$, we
find that $R\approx F$ (see main text and Fig. 16). Moreover, it holds that
$\|R-F\|_{2}\to 0$ as $P_{ij}\to 1$ element-wise. This implies that functional
relationships between the phases $\boldsymbol{\theta}^{*}(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.

###  1.9 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 $\mathcal{V}_{1}$ and load buses $\mathcal{V}_{2}$. A structure-
preserving power network model contains $|\mathcal{V}_{1}|$ second-order
Newtonian and $|\mathcal{V}_{2}|$ first-order kinematic phase oscillators
obeying [64]:

$\displaystyle\begin{cases}M_{i}\ddot{\theta}_{i}+D_{i}\dot{\theta}_{i}&=\omega_{i}+\sum_{j=1}^{|\mathcal{V}_{1}|}a_{ij}\sin(\theta_{j}-\theta_{i}),\
i\in\mathcal{V}_{1},\\\\[3.00003pt]
D_{i}\dot{\theta}_{i}&=\omega_{i}+\sum_{j=1}^{|\mathcal{V}_{2}|}a_{ij}\sin(\theta_{j}-\theta_{i}),\
i\in\mathcal{V}_{2},\end{cases}$ (26)

where $M_{i}$, $D_{i}$ are the generator inertia constant, and the damping
coefficient, respectively. In the equation for the generators
$\mathcal{V}_{1}$, $\omega_{i}=P_{m,i}$, which is the mechanical power input
from the prime mover, and in the equation for the load buses
$\mathcal{V}_{2}$, $\omega_{i}=P_{\ell,i}$, which denotes the real power drawn
by load $i$. Finally, the weight $a_{ij}$ equals
$a_{ij}=|v_{i}||v_{j}|\mathrm{Im}({Y}_{ij})$, 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.

We assume that thermal limit constraints are equivalent to phase cohesiveness
requirements. To be precise, we obtain a bounded power flow
$a_{ij}\sin(\theta_{j}-\theta_{i})$ for the line $(i,j)$ whenever the angular
distance $|\theta_{j}-\theta_{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 $\boldsymbol{v}$, $\boldsymbol{p}_{\ell}$