Functional control of oscillator networks

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.

Fig. 1. Network control to enforce a desired functional pattern from an abnormal or undesired one.

The left panel contains a network of n = 7 oscillators (top left panel, line thickness is proportional to the coupling strength), whose vector of natural frequencies ω has zero mean. The phase differences with respect to θ₁ (i.e., θ_i−θ₁) 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 (bottom left panel, color coded). The center left panel depicts the functional pattern R corresponding to such phase differences over time. The right panel illustrates the same oscillator network after a selection of coupling strengths and natural frequencies have been tuned (in red, the structural parameters A and ω are adjusted to A_c and ω_c) to obtain the phase differences $\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 center right. In this example, we have computed the closest set (in the ℓ₁-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.

Fig. 2. Mapping between desired phase differences and interconnection weights.

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

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 $\overline{B}$ and vector of natural frequencies ω. The right side illustrates the cone generated by the columns of $\overline{B}$. In this example, $\mathcal{S}=\left\{1,2\right\}$ satisfies the conditions for the existence of δ ≥ 0 in Eq. (5), as ω is contained within the cone generated by the columns ${\overline{B}}_{:,\mathcal{S}}$. b The (directed) Hamiltonian path described by the columns of ${\overline{B}}_{:,\mathcal{H}}$, with $\mathcal{H}=\left\{1,2\right\}$, in the network of panel (a). c The existence of such an Hamiltonian path, together with a positive projection of ω onto ${\overline{B}}_{:,\mathcal{H}}$, also ensure that there exists a strictly positive δ > 0 solution to BD(x)δ = ω. In particular, for any choice of x₁₂, x₂₃ ∈ (0, π), Eq. (4) reveals that if $0<A_{13}<0.5/\sin \left(x_{12}+x_{23}\right)$, then there exist strictly positive weights A₁₂ > 0 and A₂₃ > 0 such that δ > 0.

Fig. 4. The intersection of an affine space with sin(xdep) determines the compatible functional patterns of 3 identical oscillators.

Consider a fully connected network of n = 3 identical oscillators with zero natural frequency and δ = 1. It is well known that x⁽¹⁾ = [0 0]^T, x⁽²⁾ = [π 0]^T, x⁽³⁾ = [0 π]^T, x⁽⁴⁾ = [π π]^T are phase difference equilibria. Furthermore, because $\sin \left(\theta \right)=\sin \left(\pi -\theta \right)$, this figure illustrates $\sin \left(x_{13}\right)$ as a function of x₁₂ and x₂₃ in four different panels: $\sin \left(x_{12}+x_{23}\right)$ (top left), $\sin \left(\pi -x_{12}+x_{23}\right)$ (top right), $\sin \left(x_{12}+\pi -x_{23}\right)$ (bottom left), and $\sin \left(-x_{12}-x_{23}\right)$ (bottom right). The fourth panel reveals that the two functional patterns compatible with x^(j), j ∈ {1, …, 4}, correspond to x⁽⁵⁾ = [2π/3 2π/3]^T and x⁽⁶⁾ = [−2π/3 −2π/3]^T (in red).

Fig. 5. A homogeneous cycle network admits infinite compatible functional patterns.

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

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 ${\mathit{x}}_{\mathrm{desired}}={\left[\frac{4\pi }{5}\frac{\pi }{3}\frac{\pi }{10}\right]}^{\mathsf{T}}$. The starting pattern ${\mathit{x}}_{\mathrm{original}}={\left[\frac{\pi }{10}\frac{\pi }{3}\frac{4\pi }{5}\right]}^{\mathsf{T}}$ is associated with interconnection weights δ = [3.4026 3.4641 6.4721]^T. Applying the optimal correction α^* yields positive interconnection weights δ + α^* = [6.4721 3.4026 3.4641]^T that achieve the desired functional patterns x_desired. b Joint allocation of two compatible equilibria for the phase difference dynamics. By taking θ₁ as a reference, we choose two points for the phase differences x_1i = θ_i − θ₁, i ∈ {2, …, 7}, to be imposed as equilibria in a network of n = 7 oscillators: ${\mathit{x}}_{\mathrm{desired}}^{\left(1\right)}={\left[\frac{\pi }{6}\frac{\pi }{6}\frac{\pi }{4}\frac{\pi }{4}\frac{\pi }{6}\frac{\pi }{4}\right]}^{\mathsf{T}}$ and ${\mathit{x}}_{\mathrm{desired}}^{\left(2\right)}={\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 (δ + α^*) that solves the minimization problem (11) with constraint (12). The trajectories start at the (unstable) equilibrium point ${\mathit{x}}_{\mathrm{desired}}^{\left(1\right)}$ at time t = 0, and converge to the point ${\mathit{x}}_{\mathrm{desired}}^{\left(2\right)}$ after a small perturbation $\mathit{p}\in {\mathbb{T}}^7$, with π ∈ [0 0.05], is applied to the phase difference trajectories at time t = 50.

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 (13) to achieve stability of the pattern ${\mathit{x}}_{\mathrm{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 (λ₂ = 0.0565, in red). The optimal correction α^* is gradually applied to the existing interconnections from the left-most panel to the right-most one at $\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 (λ₁ = 0 due to rotational symmetry and λ₂ = − 0.0178, in green).

Fig. 8. Replication of empirically recorded functional connectivity in the brain through tuning of the natural frequencies of Kuramoto oscillators.

The anatomical brain organization provides the network backbone over which the oscillators evolve. The filtered fMRI time series provide the relative phase differences between co-fluctuating brain regions, and thus define the desired phase differences x, which is used to calculate the metabolic change encoded in the oscillators' natural frequencies. In this figure, we select the 40-s time window from t₀ = 498 s to t_f = 538 s for subject 18 in the second scanning session. We obtain ∥R−F∥₂ = 0.2879. Additionally, we verify that the analysis of the Jacobian spectrum (see Eq. (10)) accurately predicts the stability of the phase-locked trajectories. Supplementary Fig. 7a illustrates the basin of attraction of R, which we numerically estimate to be half of the torus.

Fig. 9. Fault recovery in the IEEE 39 New England power distribution network through minimal and local intervention.

a New England power distribution network. The generator terminal buses illustrate the type of generator (coal, nuclear, hydroelectric). We simulate a fault by disconnecting the transmission line 25 (between loads 13 and 14). b The fault causes the voltage phases θ to depart from normal operating conditions, which could cause overheating of some transmission lines (due to violation of the nominal thermal constraint limits) or abnormal power delivery. To recover the pre-fault active power flow and promote a local (sparse) intervention, we solve the optimization in Eq. (11) by minimizing the 1-norm of the structural parameter modification δ with no scaling parameters in the cost functional. The network returns to the initial operative conditions with a localized modification of the neighboring transmission lines' impedances.