Functional Control of Oscillator 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.

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}$.

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.

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.

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.

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.

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.

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).

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.

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

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

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.

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.

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}$.

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.

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}}}$.

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}$.

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.

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.

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

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.

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.

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}$: