Emergent Topology of Optimal Networks for Synchrony

TL;DR

A gradient-based optimization framework is developed to identify synchrony-optimal weighted networks under a constrained coupling budget, which exhibit counterintuitive features: they are sparse, bipartite, elongated, and extremely monophilic.

Abstract

Real-world networks, whether shaped by evolution or intelligent design, are typically optimized for functionality while adhering to physical, geometric, or budget constraints. Yet tools to identify such structures remain limited. We develop a gradient-based optimization framework to identify synchrony-optimal weighted networks under a constrained coupling budget. The resulting networks exhibit counterintuitive features: they are sparse, bipartite, elongated, and extremely monophilic (i.e., the neighbors of any node are similar to one another while differing from the node itself). These findings are matched by "constructive" theory: a nonlinear differential equation identifies which pairs of nodes are coupled, while a variational principle prescribes the budget allocated to each node. Dynamics unfolding over optimal networks provably lack a synchronization threshold; instead, as the budget exceeds a calculable critical value, the system globally phase-locks, exhibiting critical scaling at the transition. Our results have implications for power grids, neuromorphic computing, and other coupled oscillator technologies.

Figures (5)

• Figure 1: Optimized network for Kuramoto synchronization under a limited coupling budget $b$. A representative network is shown on top, while snapshots at successive stages of optimization are visualized below. Node colors encode intrinsic frequencies $\omega_i$; edge thicknesses reflect edge weights $A_{ij}>0$. The emergent structure is highly sparse and bipartite, with edges exclusively linking nodes of contrasting frequencies (e.g., blue to red). It is surprisingly elongated, with long average path lengths. The network also exhibits extremely strong monophily: the neighbors of any node are strikingly homogeneous. This effect is readily observed at the periphery. For example, on the left, dark-red nodes have neighbors that are all shaded light-blue.
• Figure 2: Computational graph for gradient-based network optimization.Left: Forward pass. A tunable parameter matrix $\boldsymbol{P}$ defines a symmetric adjacency matrix $\boldsymbol{A}$ obeying the budget constraint. The dynamical system \ref{['oscillator_eq']} is then numerically integrated from an initial condition $\boldsymbol{\theta}_0=[\theta_1(0),\dots,\theta_N(0)]^T$, where each subsequent state $\boldsymbol{\theta}_t$ depends on the previous state $\boldsymbol{\theta}_{t-1}$ and on $\boldsymbol A$. We compute $r_t$ at each time step and then the objective function $\langle r \rangle$ over a duration $\mathcal{T}$. An initial time window may be discarded to ignore the initial transient behavior. Right: Backward pass. The dependency structure of $d\langle r \rangle / dP_{ij}$ is illustrated via so-called "backpropagation"---that is, the chain rule of calculus. For clarity, we use color to highlight only the backward pathways passing through $r_4$, which correspond explicitly to the term $\frac{\partial \langle r \rangle}{\partial r_4}\frac{d r_4}{d P_{ij}}$. The computational graph and chain-rule are automatically tracked and computed using standard computational techniques for machine learning paszke2017automatictorchdiffeq. The code is openly available on GitLabmikaberidze2025networkoptimization.
• Figure 3: Disappearance of the synchronization threshold and emergence of global phase-locking.(Left) Time series of phases $\theta_i(t)$ for a synchrony-optimized Kuramoto network (and final phase snapshots on the unit circle) are shown for a varying coupling budget $b$. A cluster of phases exists for all $b$, and global phase-locking occurs above a critical budget $b_c$. (Top-right) We plot the time-averaged Kuramoto order parameter and size of the largest phase-locked cluster versus $b$ for both a fixed all-to-all network (orange) and optimal networks (blue). Optimal networks lack a synchronization threshold, and $b_c$ is shown by the dashed black line given by Eq. \ref{['b_c']}. (Bottom-right) Phase diagram showing the absence/presence of global phase-locking for optimal networks as a function of coupling budget $b$ and a frequency distribution parameter $\alpha$. Color indicates the outcome of numerical optimization and simulation. In all panels, $N=100$ intrinsic frequencies $\omega \in [-1, 1]$ are drawn from the family $g_\alpha(\omega) \propto 1/(1 + \alpha \omega^2)$. For this family, Eq. \ref{['b_c']} yields $b_c = \log(1+\alpha) / (2\sqrt{\alpha} \arctan\sqrt{\alpha})$. The bottom-right diagram spans $\alpha \in [-1, 1]$, interpolating between bimodal and unimodal frequency distributions, while the other panels use fixed $\alpha = 1$.
• Figure 4: Edge weights $A_{ij}$ as a function of oscillator frequencies$\omega_i$ and $\omega_j$ in a synchrony-optimized Kuramoto network. Black curves in either subplot mark the negative solution $\omega_j = \nu_-(\omega_i)$ of Eq. \ref{['curve_eq']}. (a) Each blue curve traces the distribution of neighbor frequencies for a given $\omega_i$, with density peaking near the predicted $\nu_-$. Results are shown after 2000 optimization steps. (b) A heatmap representation of the network after 5000 additional steps, where the alignment sharpens significantly within the frequency$\times$frequency space and closely matches our analytical prediction $\nu_-(\omega)$ (dashed black curves). Both panels use $N=10^4$, $b=1$, and $g_1(\omega) \propto 1/(1 + \omega^2)$.
• Figure 5: Analytical and numerical results for the strong-coupling regime.(a) Kuramoto order parameter $r$ from Eq.\ref{['eq:strongly_coupled:r']} evaluated using the two pairing function branches $\nu_+(\omega)$ and $\nu_-(\omega)$ across the family of frequency distributions $g_\alpha(\omega) \propto 1/(1 + \alpha \omega^2)$. The nontrivial branch $\nu_-$ consistently outperforms $\nu_+$, explaining the convergence of our gradient-based optimization to $\nu_-$. (b) Node strengths $s(\omega)$ from Eq.\ref{['eq:strongly_coupled:s']} (dashed) compared with optimized results (solid) after $10^5$ training steps with $N=10^4$, $b=20$, and $\alpha=1$. (c) Stationary phases $\theta^*(\omega)$ from Eq. \ref{['eq:strongly_coupled:theta']} (dashed) versus simulation on the numerically optimized network (solid). These reported $r$ values indicate near-perfect agreement between analytical and numerical solutions, with the analytical prediction being slightly higher than the numerics.