Random heterogeneity outperforms design in network synchronization

We consider a network of $N$ delay-coupled nonidentical Stuart–Landau oscillators, whose dynamics are governed bywhere $z_j=r_j{e}^{\mathrm{i}{\psi }_j}=x_j+\mathrm{i}{y}_j$ is a complex variable representing the state of the $j$th oscillator and $K$ is the coupling strength. The adjacency matrix $\mathit{A}=\left\{A_{jk}\right\}$ encodes the network structure, and $d_j={\sum }_k{A}_{jk}$ is the indegree for oscillator $j$. The coupling delay $\tau$ models the finite speed of signal propagation in real systems, which is often significant in biological (28–30), physical (31–33), and control systems (34, 35). The local dynamics $f_j$ are expressed in the following canonical form for systems born out of a Hopf bifurcation (36):where ${\lambda }_j$, ${\omega }_j$, and ${\gamma }_j$ are real parameters associated with the amplitude, base frequency, and amplitude-dependent frequency of the underlying limit-cycle oscillations.

The oscillators are identical when ${\lambda }_j=\lambda ,\hspace{0.17em}{\omega }_j=\omega$, and ${\gamma }_j=\gamma$ for all $j$. For identical oscillators, the identical synchronization state is defined by the limit-cycle synchronous solutionThe amplitude $r_0$ and angular velocity ${\mathrm{\Omega }}_0$ can be found by solving the transcendental equationswhere $\mathrm{\Phi }=-{\mathrm{\Omega }}_0\tau$ (10). When random heterogeneity is introduced through one or more oscillator parameters, the identical synchronization state [3] will, in general, no longer exist. Nonetheless, we show that heterogeneous systems can still admit states that are synchronized in the sense of exhibiting cohesive phase and amplitude dynamics, as formalized below. Here, we consider synchronization in this sense and ask whether it can be stabilized by random oscillator heterogeneity.

We start by considering a homogeneous system consisting of $N=18$ identical Stuart–Landau oscillators coupled through a directed ring network ($A_{jk}=1$ if $k=j+1\mathrm{mod}N$ and $A_{jk}=0$ otherwise) for $\lambda =0.1$, $\omega =-0.28$, $\gamma =-4.42$, $K=0.3$, and $\tau =1.8\pi$. Under this parameter choice, the limit-cycle synchronous solution [3] is unstable, and the system evolves into a symmetry-broken state, exhibiting incoherent chaotic dynamics. As an example of a heterogeneous system, we consider the same network with the base angular velocity ${\omega }_j$ drawn from a Gaussian distribution with mean $\omega =-0.28$ (as in the homogeneous system) and SD $\sigma =0.1$. For additional details on the numerical procedure, see Materials and Methods.

In Fig. 1, we show typical trajectories, order parameters, and angular velocities for the homogeneous and heterogeneous systems. Here, two order parameters are introduced to measure the cohesiveness of the dynamics: the phase order parameter $R_1\hspace{0.17em}=\hspace{0.17em}\left|{\sum }_j{e}^{\mathrm{i}{\psi }_j}/N\right|$, which is the one typically used in the study of Kuramoto oscillators; and the phase-amplitude order parameter $R_2=\frac{1}{\text{max}\left(r_j\right)}\left|{\sum }_j{r}_j\hspace{0.17em}{e}^{\mathrm{i}{\psi }_j}/N\right|$, which measures the coherence in both phases and amplitudes. It follows that both $R_1$ and $R_2$ are constant for frequency-synchronized states. For the trajectories shown, the two systems were initialized close to the limit-cycle synchronous state. The homogeneous system loses synchrony at t = 1,500 and transitions to an incoherent state with both $R_1$ and $R_2$ fluctuating around 0.2. Remarkably, despite having different base frequencies, oscillators in the heterogeneous system converge to a stable cohesive state with large order parameters ($R_1>R_2>0.9$) and identical angular velocities. That is, the oscillators are not only approximately synchronized in phase and amplitude—they are also exactly synchronized in frequency (i.e., phase-locked). For an animation of this phenomenon in larger systems, see Movie S1. It is worth noting that the same phenomenon is also observed when the coupling function is nonlinear and/or when the coupling delay is link-dependent, which we demonstrate in SI Appendix, SI Text and Figs. S1 and S2.

To gain theoretical understanding of the effect shown in Fig. 1 and its prevalence, we characterize the synchronization states of interest in the presence of heterogeneity and derive analytical conditions for their stability. These results are established for delay-coupled Stuart–Landau oscillators with arbitrary heterogeneity.

Inspired by the fact that the heterogeneous system in Fig. 1B settles into a frequency-synchronized state, for which the frequencies of all oscillators are equal and the phase differences and amplitudes remain constant, we employ the following ansatz:where oscillator $j$ has amplitude $r_j$ and phase lag ${\delta }_j$, both of which are assumed to be constant, and all oscillators share the same angular velocity $\mathrm{\Omega }$. Substituting this ansatz into Eq. 1, we obtain $2N$ nonlinear algebraic equations with $2N$ unknowns:for $j=1,\dots ,N$, where ${\mathrm{\Phi }}_{jk}={\delta }_k-{\delta }_j-\mathrm{\Omega }\tau$. Taking ${\delta }_1=0$, which can be done without loss of generality, the solution of Eqs. 6a and 6b determines ${\delta }_2,\dots ,{\delta }_N,r_1,\dots ,r_N$, and $\mathrm{\Omega }$. As shown in Fig. 2, when parameter heterogeneity is not too large, this gives us frequency-synchronized states that are close to the identical synchronization state of the homogeneous system given by Eqs. 3 and 4.

We can analyze the stability of the frequency-synchronized states through the variational equation that governs the evolution of small deviations $\delta r_j\left(t\right)$ and $\delta {\psi }_j\left(t\right)$ from those states. Taking $z_j\left(t\right)=r_j\left[1+\delta r_j\left(t\right)\right]e^{\mathrm{i}\left[\mathrm{\Omega }t+{\delta }_j+\delta {\psi }_j\left(t\right)\right]}$, ${\mathit{\eta }}_j={\left(\delta r_j,\delta {\psi }_j\right)}^{⊺}$, and $\mathit{\eta }={\left({\mathit{\eta }}_1^{⊺},\dots ,{\mathit{\eta }}_N^{⊺}\right)}^{⊺}$, the variational equation can be written aswhere ${\mathit{J}}_j=\left(\begin{array}{cc}-2r_j^2& 0\\ -2{\gamma }_j{r}_j^2& 0\end{array}\right)$, $\mathit{D}=\text{diag}\left(\frac{1}{d_1},\dots ,\frac{1}{d_N}\right)$, ${\mathit{R}}_{jk}=\frac{r_k}{r_j}\left(\begin{array}{cc}\cos {\mathrm{\Phi }}_{jk}& -\sin {\mathrm{\Phi }}_{jk}\\ \sin {\mathrm{\Phi }}_{jk}& \cos {\mathrm{\Phi }}_{jk}\end{array}\right)$, ${\mathit{P}}_j=\frac{1}{d_j}{\sum }_k{A}_{jk}{\mathit{R}}_{jk}$, and $\oplus \left({\mathit{J}}_j-K{\mathit{P}}_j\right)=\text{diag}\left({\mathit{J}}_1-K{\mathit{P}}_1,\dots ,{\mathit{J}}_N-K{\mathit{P}}_N\right)$. Here, $\mathit{D}\mathit{A}\otimes {\mathit{R}}_{jk}$ is a $2N\times 2N$ matrix obtained by replacing each entry $A_{jk}/d_j$ in $\mathit{D}\mathit{A}$ with a $2\times 2$ block given by $A_{jk}{\mathit{R}}_{jk}/d_j$.

Since all of the matrices in Eq. 7 are time-independent, we can assume $\mathit{\eta }\left(t\right)=\mathit{\eta }\left(0\right)e^{v_{\ell }t}$ and reduce the stability calculation to determining the exponents $v_{\ell }$ that solve the following characteristic equation:where $\ell$ indexes the solutions. The Lyapunov exponents of Eq. 7 are given by the real parts of $v_{\ell }$. One Lyapunov exponent (referred to as $\text{Re}\left(v_0\right)$) is identically null and corresponds to the perturbation mode parallel to the frequency-synchronization manifold (in which the phases of all oscillators are subject to the same perturbation). The other Lyapunov exponents correspond to perturbation modes transverse to the synchronization manifold. The stability of the frequency-synchronization state is determined by the sign of the maximum transverse Lyapunov exponent (MTLE), which is given by $\mathrm{\Lambda }=\underset{\ell \ne 0}{max}\left\{\text{Re}\left(v_{\ell }\right)\right\}$.

We now examine systematically the phenomenon of synchronization induced by random heterogeneity. In particular, we address key questions underlying its prevalence. For example, what is the effect of the magnitude of parameter mismatches? Do the results change significantly depending on which parameters are made heterogeneous? Most importantly, can different realizations of random heterogeneity consistently induce synchronization?

In Fig. 3, we start with the same homogeneous system as in Fig. 1A and introduce heterogeneity in $\left\{{\omega }_j\right\}$, $\left\{{\lambda }_j\right\}$, and $r_0^2\left\{{\gamma }_j\right\}$, respectively. (Here, the factor $r_0^2$ is introduced to scale $\sigma$ for ${\gamma }_j$, because the influence of ${\gamma }_j$ in Eq. 2 is scaled by the square of the oscillation amplitude. The constant $r_0$ can be found by solving Eqs. 4a and 4b for the corresponding homogeneous system.) In all cases, the SD is $\sigma$, and the mean is taken to be the same as the corresponding parameter in the homogeneous system. For each realization of heterogeneity, as $\sigma$ increases from zero, the identical synchronization state progressively changes into a phase-locked state with large order parameters. The stability of this state is measured by $\mathrm{\Lambda }\left(\sigma \right)$, which we obtain by solving Eq. 8 for each realization of heterogeneity. The filled green curves in Fig. 3 A–C, Upper, show the probability that synchronization is stabilized by random heterogeneity in each parameter. These results are verified by direct simulations of Eqs. 1 and 2 for various $\sigma$, shown as purple circles. In Fig. 3 A–C, Lower, we plot $\mathrm{\Lambda }\left(\sigma \right)$ for a representative subset of realizations of heterogeneity in each parameter, visualizing their impact on stability as an ensemble.

One can see from Fig. 3 that there is always a sweet spot of optimal heterogeneity at an intermediate value of $\sigma$. Around that sweet spot, the green curves stay very close to one, indicating that intermediate levels of heterogeneity can consistently induce synchronization, largely independent of its particular realization. It is interesting to note from Fig. 3 A–C, Lower, that small heterogeneity always improves stability under the conditions considered, as reflected in the monotonic decrease of $\mathrm{\Lambda }\left(\sigma \right)$ for small $\sigma$. Disorder can also consistently stabilize synchronization when all three parameters are allowed to be heterogeneous, as demonstrated in SI Appendix, SI Text and Fig. S3. Furthermore, we verified that the same effect can be observed for a wide range of network sizes and different network structures (SI Appendix, SI Text and Figs. S4–S6).

It is important to compare the effect of random and nonrandom heterogeneities. When the heterogeneity is purposively designed, Stuart–Landau oscillators can synchronize identically (i.e., all phase differences are identically zero and all amplitudes are equal), even though they are nonidentical. This is most easily seen from Eqs. 4a and 4b, whose solution remains invariant under the transformation $\omega \to \omega +h$, $\gamma \to \gamma +h/r_0^2$ for any $h\in \mathbb{R}$. Thus, any given Stuart–Landau oscillator belongs to a continuous family of nonidentical Stuart–Landau oscillators parameterized by $h$, within which the oscillators can synchronize identically with each other. Moreover, as shown in ref. 25, mixing different oscillators from the same family can stabilize identical synchronization that would otherwise be unstable.

By designing heterogeneity to preserve identical synchronization, can we do better than by relying on random heterogeneity? Once again, we start with the homogeneous system studied in Figs. 1 and 3. The oscillators are then made heterogeneous by sampling from the identically synchronizable family, with $h$ drawn from a Gaussian distribution. More concretely, ${\omega }_j=\omega +h_j$ and ${\gamma }_j=\gamma +h_j/r_0^2$, where $\left\{h_j\right\}$ has SD $\sigma$ and mean zero. This can be seen as a special subset of oscillators with random heterogeneity in both parameters $\left\{{\omega }_j\right\}$ and $\left\{{\gamma }_j\right\}$, the crucial difference being that ${\omega }_j-\omega$ and ${\gamma }_j-\gamma$ are not independent for the designed heterogeneity.

In Fig. 4, we compare the ensemble average MTLE and order parameters between systems with random heterogeneity and systems with designed heterogeneity on the same parameters. Consistent with Fig. 3, random heterogeneity is most effective for intermediate magnitudes $\sigma$, ranging from 0.05 to 0.1. On the other hand, designed heterogeneity is effective for much larger $\sigma$, from about 0.4 to 0.6, which may be interpreted as a consequence of the identical synchronization solution being preserved in this case. Remarkably, no system with designed heterogeneity is stable within the range for which random heterogeneity is effective. This implies that at intermediate magnitude, random heterogeneity can outperform heterogeneity specifically designed to preserve identical synchronization.

To gain further understanding, in Fig. 5, we focus on a minimal system formed by three nonidentical Stuart–Landau oscillators coupled through a directed ring network. The $j$th oscillator has parameters $\left\{{\lambda }_j,{\omega }_j,{\gamma }_j\right\}=\left\{\lambda ,\omega +h,\gamma +\left(h+{\mathrm{\Delta }}_j\right)/r_0^2\right\}$, with the constraint that ${\sum }_{j=1}^3{\mathrm{\Delta }}_j=0$. The parameter $h$ is introduced to vary the synchronization stability of the homogeneous system without altering the synchronous solution. This enables us to investigate all possible realizations of heterogeneous ${\gamma }_j$ for different levels of instability by sweeping the ${\mathrm{\Delta }}_1$–${\mathrm{\Delta }}_2$ plane.

In Fig. 5 A and B, the origin is the only point corresponding to a homogeneous system, and the differences among oscillators increase as one moves away from the origin along the radial directions. Stability analysis indicates that regions of stability appear for intermediate levels of oscillator heterogeneity, as shown in Fig. 5A ($\mathrm{\Lambda }<0$, blue belts). The phase-locked state is unstable for weak disorder ($\mathrm{\Lambda }>0$, red areas) and ceases to exist for strong disorder (blank areas). A complementary perspective is offered by direct simulations, as shown in Fig. 5B. Because order parameters averaged over time are a poor indicator of coherence for systems with a small number of oscillators, we quantify the level of coherence using the minimum of $R_2$ over a period of 10,000 time units after the initial transient. For zero and small heterogeneity, the three oscillators are in an incoherent state, with $min{R}_2\approx 0$. As $\sigma$ is increased further, the oscillators first settle into an approximate synchronization state with $min{R}_2$ ranging from 0.6 to 0.9 (light purple regions). The level of coherence continues to improve until it plateaus at $min{R}_2\approx 0.96$ for phase-locked states (dark purple regions), which correspond to the stable states marked by the blue belts in Fig. 5A. Finally, once we cross the outer boundary, synchrony is lost again, and the value of $min{R}_2$ falls back to approximately zero. This incoherence–coherence–incoherence transition is illustrated in Fig. 5C with representative trajectories from each stage. It is worth noting that even before the phase-locked state is fully stabilized, disorder can already induce approximate synchronization states with well-defined rhythms, as illustrated by the second trajectory.

Fig. 5 demonstrates two competing effects of disorder: When heterogeneity is too small, it cannot tame synchronization instability; when it is too large, it destroys the synchronization state. In other words, there is a trade-off between synchronizability and stability, and stable synchronization naturally emerges at intermediate levels of disorder. Another interesting observation is that the stable belts are contiguous in all cases in Fig. 5A and completely surround the unstable regions in the middle, which explains why intermediate levels of heterogeneity can consistently stabilize synchronization. It also demonstrates that the effect is robust against increasing instability (controlled by $h$) in the homogeneous system.

A natural question at this point is whether the described phenomenon is robust and general enough to be observed in real systems. To provide an answer, we performed experiments using chemical oscillators based on the electrochemical dissolution of nickel in sulfuric acidic medium (26). The experimental apparatus consists of a counter electrode, a reference electrode, a potentiostat, and $N$ nickel wires submerged in the same sulfuric acidic medium, each attached to a resistor (Fig. 6A). At constant circuit potential ($V_0=1.24$ V relative to the reference electrode) and with the resistance of resistors set to ${\xi }_j=1.06$ kohm, the dissolution rate of each nickel wire, measured as its current, exhibits periodic oscillations (37). The oscillatory dynamics originate from a Hopf bifurcation at $V_0=1.07$ V. Compared to the circuit potential used in some previous studies (33), the system here is farther away from the bifurcation point. When the wires are placed sufficiently far from each other, the current oscillations do not show noticeable synchronization, confirming that the interactions through the solution are negligible. Coupling among the wires can be introduced through external feedback (26, 33), in which the circuit potentials of the wires $V_j\left(t\right)$ are set based on the measured currents $I_j\left(t\right)$ asfor $j=1,\dots ,N$, where $K$ and $\tau$ are the experimental coupling strength and delay, respectively. Here, we investigate $N=16$ wires with oscillatory currents arranged in an undirected four-by-four lattice network with periodic boundary conditions, which can be seen as a two-dimensional variant of the ring networks considered above. For additional details on the experimental setup and procedure, see Materials and Methods.

For relatively strong coupling ($K\approx -0.40$ V/mA) and no delay ($\tau =0$ s), the system exhibits in-phase synchronization (38). Similar in-phase synchronization exists for large delay ($\tau \approx 2.4$ s) that corresponds to the mean period of the uncoupled oscillations. When $\tau$ is set to $\approx 1.2$ s (about half of the oscillation period), the system exhibits a two-cluster state in which every other element on the grid is in phase, and the neighboring elements are in antiphase. When the delay is set between these two regions ($\tau \approx 1.75$ s), the system exhibits a desynchronized state. Nominal oscillator heterogeneity was introduced by setting the resistance of each oscillator to a different value ${\xi }_j$, while keeping the mean resistance fixed to $\overline{\xi }=1.06$ kohm. The level of nominal heterogeneity is measured by the SD $\sigma$ among all ${\xi }_j$.

First, we randomly picked one realization of heterogeneity and experimentally tested its effect on the collective dynamics at different levels of heterogeneity $\sigma$. Each experiment was initiated close to the in-phase synchronization state and consisted of running the system for 600 s. The level of coherence was measured by the synchronization error e$\left(t\right)$, defined as the SD among the currents $I_j$ at time $t$. In this case, the synchronization error is a more natural measure of coherence than order parameters because the experimental system is not in the close vicinity of a Hopf bifurcation, and the dynamics of the amplitude variables are oscillatory. (Nevertheless, we verified that the order parameters of the phases extracted using either Hilbert transform or peak-detection algorithms give similar results as the ones obtained using the synchronization error.) The experimental results summarized in Fig. 6B reveal a well-defined minimum of the average synchronization error ⟨e⟩ (averaged over the last 200 s of each experimental run) for an intermediate level of nominal heterogeneity, $\sigma =0.1$ kohm. This optimization of synchronization at intermediate levels of heterogeneity is consistent with what we observed numerically for delay-coupled Stuart–Landau oscillators.

Unlike the idealized systems used in simulations, experimental systems come with unavoidable imperfections and uncertainties. As a result, the electrochemical oscillators in our experiments have slightly different dynamics, even when the resistances are all set to the same nominal value. These relatively small inherent heterogeneities can arise because of unavoidable differences in the metal wires (e.g., in composition and size) and surface conditions (oxide film layer thickness, localized corrosion, etc.). To account for such inherent heterogeneity, we use peak-detection algorithms (39) to extract the natural frequency and amplitude of each uncoupled oscillator, and we use that information to calculate the measured oscillator heterogeneity $\mathrm{\Delta }$ for both systems with homogeneous ${\xi }_j$ and systems with heterogeneous ${\xi }_j$. Here, $\mathrm{\Delta }={\overline{\rho }}_T+{\overline{\rho }}_A$, where ${\overline{\rho }}_T$ (${\overline{\rho }}_A$) is the SD of the oscillation periods (amplitudes) of the uncoupled oscillators normalized by the mean. (For additional details on the data-analysis protocol, see Materials and Methods.) In Fig. 6C, we show results for five sets of independent experiments. Each experiment corresponds to a different realization of heterogeneous resistances (for $\sigma$ fixed at 0.13 kohm) and of the homogeneous system (corresponding to $\sigma =0$ kohm). It can be seen that, when uncoupled, all heterogeneous systems have a much higher measured oscillator heterogeneity $\mathrm{\Delta }$ than the homogeneous systems. In contrast, when coupled, the heterogeneous systems achieve significantly better coherence than the homogeneous systems, which is reflected by a consistently smaller ⟨e⟩.

The striking difference between the behavior of the homogeneous and heterogeneous systems is further visualized in Fig. 7. There, we compare the dynamics in the first ($\sigma =0$ kohm) and the fifth ($\sigma =0.13$ kohm) data points from Fig. 6B. The time series of the homogeneous system (Fig. 7A) is very much incoherent compared to that of the heterogeneous system (Fig. 7B). Accordingly, the heterogeneous system exhibits a smaller synchronization error and a more regular rhythm, as shown in Fig. 7 C and D. In SI Appendix, SI Text and Fig. S7, we show that the same phenomenon can be observed when random shortcuts are added to the four-by-four lattice [which creates a small-world network (40)], suggesting that our conclusions extend to networks with heterogeneous degrees.

Delay-coupled Stuart–Landau oscillators were simulated by employing the dde23 integrator in MATLAB, with the relative and absolute tolerances both set to $10^{-4}$. To initialize an oscillator network, we introduced a random perturbation of the order of $10^{-1}$ to the synchronization state at $t=0$. Each system was then evolved for 10,000 time units, which is long enough for the oscillators to settle into either a coherent state (if synchronization is stable) or an incoherent state (if synchronization is unstable). Our code for simulating delay-coupled Stuart–Landau oscillators can be found at https://github.com/y-z-zhang/disorder_sync.

The experiments were performed by using a standard three-electrode cell with a platinum counter, a Hg/Hg₂SO₄/sat$\cdot$K₂SO₄ reference, and a nickel-array working electrode. The electrolyte was 3 M H₂SO₄ at 10 °C. The electrode array consisted of 16 1-mm-diameter nickel wires with a spacing of 3 mm. The wires were embedded in epoxy, so that only the wire ends were exposed to the electrolyte. Before the experiments, the electrode array was polished with a series of sandpapers. A multichannel potentiostat (catalog no. Gill-IK64, ACM Instruments), interfaced with a real-time LabVIEW controller (60), was used to measure the current $I_j\left(t\right)$ and set the potential $V_j\left(t\right)$ of the $j$th wire according to Eq. 9 at a rate of 200 Hz. Throughout the experiments, we set the circuit potential to $V_0=$ 1.24 V. Without heterogeneity, the individual resistors were set to 1.06 kohm. Heterogeneity was introduced by setting the individual resistors to different nominal values drawn from a normal distribution, while keeping the mean resistance fixed at 1.06 kohm. To avoid accidentally balancing out the inherent heterogeneity, only random realizations of nominal heterogeneity that had a negligible correlation with the natural frequencies of the unperturbed oscillators were used (we required that the absolute value of the correlation coefficient be smaller than 0.2). The coupling delay $\tau$ was set to 75% of the mean natural period of the oscillations, which corresponds to $\tau$ in the range of 1.50 to 1.75 s throughout the experiments. The coupling strength $K$ was set to values about 10% larger in magnitude than the desynchronization threshold (between $-0.48$ and $-0.40$ V/mA in the reported experiments).

The peak-detection algorithm finds all local maxima by comparing the neighboring values in a time series. The mean of the detected peaks is taken as the oscillation amplitude, and the mean distance between consecutive peaks is the oscillation period. Our data-analysis scripts and experimental data are available at https://github.com/y-z-zhang/disorder_sync. By following the Jupyter Notebooks included in the GitHub repository, one can explore the data interactively and reproduce the results presented in Figs. 6 and 7 and SI Appendix, Fig. S7.