Sensitive Dependence of Optimal Network Dynamics on Network Structure

The relation between network structure and dynamics is determinant for the behavior of complex systems in numerous domains. An important long-standing problem concerns the properties of the networks that optimize the dynamics with respect to a given performance measure. Here, we show that such optimization can lead to sensitive dependence of the dynamics on the structure of the network. Specifically, using diffusively coupled systems as examples, we demonstrate that the stability of a dynamical state can exhibit sensitivity to unweighted structural perturbations (i.e., link removals and node additions) for undirected optimal networks and to weighted perturbations (i.e., small changes in link weights) for directed optimal networks. As mechanisms underlying this sensitivity, we identify discontinuous transitions occurring in the complement of undirected optimal networks and the prevalence of eigenvector degeneracy in directed optimal networks. These findings establish a unified characterization of networks optimized for dynamical stability, which we illustrate using Turing instability in activator-inhibitor systems, synchronization in power-grid networks, network diffusion, and several other network processes. Our results suggest that the network structure of a complex system operating near an optimum can potentially be fine-tuned for a significantly enhanced stability compared to what one might expect from simple extrapolation. On the other hand, they also suggest constraints on how close to the optimum the system can be in practice. Finally, the results have potential implications for biophysical networks, which have evolved under the competing pressures of optimizing fitness while remaining robust against perturbations.

Popular Summary

Complex systems with many interacting components can be represented as a network, where the components and their interactions are described graphically with nodes and links, respectively. Such systems include chemical-reaction pathways, social hierarchies, and electrical circuits. Networks are often tailored to optimize the performance of the system as a whole. Examples include neuronal networks, which may have evolved to maximize information-transfer efficiency or to minimize synaptic-wiring cost, and power grids, where both maximizing stability and minimizing cost are desirable. An important long-standing problem in the study of complex systems is thus to understand the properties of such optimized networks. We show that network optimization can lead to a sensitive dependence of the dynamics on the structure of the network.

Specifically, we demonstrate that the stability of dynamical states in optimal networks can be highly sensitive to the removal or addition of links and nodes, or to small changes in link weights. In undirected networks (where all interactions are symmetric), this sensitivity arises from discontinuous transitions caused by optimization. In directed networks, it arises from the prevalence of eigenvector degeneracy. We illustrate these findings with analyses for several different types of networks, such as networks of activator-inhibitor systems and networks of synchronizing power generators.

Our results suggest that the network structure of a complex system operating near an optimum can potentially be fine-tuned or controlled to significantly enhance stability while also setting limits on how close to the optimum the system can be in practice. These findings have potential implications for man-made networks, where optimization may inadvertently lead to fragility, and for biophysical networks, which appear to have evolved under the competing pressures of optimizing fitness while remaining robust against perturbations.

Figure 1

Directed and undirected networks optimized for the stability of the network dynamics can be sensitive to weighted and unweighted structural perturbations, respectively. The graphs schematically illustrate typical behavior for systems with sensitivity (in red boxes) and systems with no sensitivity (in blue boxes). Under weighted perturbations, we actually show that all undirected networks (including nonoptimal ones) are nonsensitive (lower blue box).

Figure 2

UCM and MCC networks with $n=18$ nodes. (a) The UCM network with $\ell =3$ groups (labeled a, b, and c) of $k=6$ nodes each. All pairs of nodes belonging to different groups are connected, while all pairs within the same group are not connected, leading to a total of $m=k^2\ell (\ell -1)/2=108$ links. (b) An MCC network constructed with the same number of nodes but with one fewer link ($m=107$) and groups of unequal sizes (labeled ${\mathrm{a}}^{\prime }$, ${\mathrm{b}}^{\prime }$, and ${\mathrm{c}}^{\prime }$, and of sizes 7, 7, and 4, respectively). Note that in this case some nodes within the same group are connected (as indicated by solid blue lines). (c) The complement of the UCM network in (a), which has $\ell k(k-1)/2=45$ links. In the complement, a node pair is connected if they are in the same group, and not connected if they are from different groups. (d) The complement of the MCC network in (b), which has 46 links. Since it has one more link than what can be accommodated by three isolated groups of size 6 [as in (c)], the minimum possible size of the largest component in the complement equals 7 in this case. Note that groups ${\mathrm{a}}^{\prime }$ and ${\mathrm{b}}^{\prime }$ have missing links (indicated by dashed blue lines), which correspond to the links within groups ${\mathrm{a}}^{\prime }$ and ${\mathrm{b}}^{\prime }$ in (b). The required increase in the size of the largest component in the complement forces ${\lambda }_2$ to decrease by one.

Figure 3

UCM and MCC networks of size $n=20$ for the maximization of smallest nonzero Laplacian eigenvalue ${\lambda }_2$. For a given link density $\varphi$, the orange dot indicates ${\lambda }_2$ for MCC networks. A UCM network, when possible for that value of $\varphi$, is indicated by a green dot. As $\varphi$ increases, the eigenvalue ${\lambda }_2$ experiences discrete jumps, corresponding to sudden changes in the structure of the network. The changes in the link configuration of the network’s complement, as well as the associated jumps in the size of their largest clusters, are illustrated in the circles. At $\varphi$ values just above and below the jumps, the complement has $m_c=M(n,k_{n,m_c})$ links and $m_c=M(n,k_{n,m_c})+1$ links, respectively.

Figure 4

Sensitive dependence of the Laplacian eigenvalue ${\lambda }_2$ on link density $\varphi$ for undirected networks of size $n=50$. Each curve indicates ${\lambda }_2$ normalized by the network size $n$, relative to the link density $\varphi$. (Plots of ${\lambda }_2/n$ itself can be found in Supplemental Material, Fig. S1 [58].) The blue curve shows the result of a single run of SA to maximize ${\lambda }_2$ with a fixed number of links. Each point on the magenta curve is the average over 1 000 realizations of the ER random networks with connection probability $\varphi$. The orange and red curves indicate the MCC networks for $n=50$ [Eq. (12)] and in the limit of $n\to \infty$ [Eq. (15)], respectively. Notice the square-root singularity on the left of points $\varphi ={\varphi }_{\ell }=(\ell -1)/\ell$, $\ell =2,3,\dots$, on the red curve. The green dot near one of these singularity points indicates the UCM network with $\ell =2$ and $k=25$, which achieves the upper bound ${\lambda }_2\le \lfloor \varphi (n-1)\rfloor$ shown by the black curve.

Figure 5

Optimization induces spectral degeneracy. The maximum $\mathrm{Re}({\lambda }_2)$ is computed exactly through symbolic calculation of eigenvalues for all directed networks of size $n=3$, 4, and 5. (a)–(c) Maximum values of $\mathrm{Re}({\lambda }_2)$ as a function of the number of directed links, $m_d$ (red dots). All the other possible values of $\mathrm{Re}({\lambda }_2)$ are indicated by blue dots. The black curves indicate the upper bound $\overline{\lambda }$ and lower bound $\lfloor \overline{\lambda }\rfloor$ for the maximum $\mathrm{Re}({\lambda }_2)$. The upper bound comes from Eq. (18), while the lower bound is derived in Sec. 3b3. (d)–(f) Probability distribution of geometric degeneracy $\beta$ for ${\lambda }_2$ of all networks of a given size (blue bars) and of the $\mathrm{Re}({\lambda }_2)$-maximizing networks (red bars). For each network size, the difference between the two distributions shows that the maximization tends to increase $\beta$.

Figure 6

Distribution of scaling exponents for the Laplacian eigenvalues of 6-node optimal directed networks under random structural perturbations. (a) Example network (a directed tree) with ${\lambda }_2=\cdots ={\lambda }_6=1$. (b) Example network with ${\lambda }_2=\cdots ={\lambda }_6=5$. Each plot shows histograms of $1/\gamma$, where $\gamma$ is the scaling exponent numerically estimated from ${\lambda }_i$ (the $i$th smallest perturbed eigenvalue) computed at 1 000 equally spaced values of $\delta$ in the intervals $[0,{10}^{-3}]$ (blue), $[0,{10}^{-4}]$ (yellow), and $[0,{10}^{-5}]$ (pink). We determine $\gamma$ by applying matlab’s built-in linear least-squares algorithm in log scale [86]. Each histogram is generated by estimating $\gamma$ for 10 000 realizations of $\mathrm{\Delta }A$, where each element of $\mathrm{\Delta }A$ (corresponding to the perturbation of the weight of an existing link or the addition of a new link with small weight) is chosen randomly from the uniform distribution on $[-1,1]$. When perturbing all the off-diagonal elements of the adjacency matrix $A$ (plots on the left-hand side of each panel), the results support $1/\gamma ={\beta }_j$ for both networks. When perturbing only the existing links (plots on the right-hand side of each panel), the scaling exponent depends on the initial network: the plots support $1/\gamma =1$ for the directed tree in (a) and $1/\gamma ={\beta }_j$ for the network in (b).

Figure 7

Optimal directed networks with various levels of sensitivity to generic weighted perturbations. Example networks $G_1$, $G_2$, $G_7$, and $G_{19}$ (each with 20 nodes and 57 links) all satisfy ${\lambda }_2=\cdots ={\lambda }_{20}=3$ (and thus are optimal) but have different geometric degeneracy. Perturbing the adjacency matrix of each network as $A+\delta \mathrm{\Delta }A$, we plot in double logarithmic scale the resulting change $|\mathrm{\Delta }\lambda (\delta )|=|\lambda (\delta )-3|$ for all 19 Laplacian eigenvalues (blue dots, many of which are overlapping). The same randomly chosen $\mathrm{\Delta }A$ is used for all four networks, where each $\mathrm{\Delta }{A}_{ij}$ is drawn uniformly from the interval $[-1,1]$ if the link exists from node $j$ to node $i$, and from [0, 1] otherwise. The perturbations thus allow small increase and decrease of the weight of existing links, as well as the addition of new links with small weight. In each network, nodes of the same color indicates the same in-degree, and a bidirectional arrow represents two directed links in opposite directions. The three nodes in the center of $G_1$ form a fully connected triangle and each of the other nodes has three in-links from these center nodes. In each plot, we observe the expected scaling behavior in Eq. (20), indicated by black lines, each labeled with the corresponding scaling exponent, $1/{\beta }_j$ ($G_1$: $\beta =1$ with 19 eigenvectors and ${\beta }_1=\cdots ={\beta }_{19}=1$; $G_2$: $\beta =2$ with 15 eigenvectors and ${\beta }_1=\cdots ={\beta }_{11}=1$, ${\beta }_{12}=\cdots ={\beta }_{15}=2$; $G_7$: $\beta =7$ with 7 eigenvectors and ${\beta }_1={\beta }_2={\beta }_3=1$, ${\beta }_4=2$, ${\beta }_5=3$, ${\beta }_6=4$, ${\beta }_7=7$; $G_{19}$: $\beta =19$ with 1 eigenvector and ${\beta }_1=19$).

Figure 8

Sensitive dependence on link density $\varphi$ in physical examples of undirected networks. (a) Exponential rate of convergence ${\mathrm{\Lambda }}_{\max }$ to a synchronous state of power-grid networks. (b) Mean finite-time convergence rate $\mu$ toward synchronization in networks of optoelectronic oscillators. (c) Critical diffusivity threshold ${ϵ}_c$ for Turing instability in networks of activator-inhibitor systems. Description of these three systems can be found in Supplemental Material [58], Secs. S1A–S1C. The orange curves indicate the values of these quantities for the MCC networks with $n=100$ constructed by the procedure described in Appendix pp3, while the red curve in (c) is the finite-$n$ approximation obtained from the asymptotic formula in Eq. (15). The blue curves indicate the corresponding values for networks found by SA. The magenta curve in (b) is the mean value for the ER random networks estimated from 1000 realizations, while in (a) and (c) the values are relative to the corresponding mean value for the ER random networks (and thus zero corresponds to the ER mean value).

Figure 9

Sensitivity to weighted perturbations in directed networks of optoelectronic oscillators. (a) Mean convergence rate $\mu$ as a function of $\delta$, illustrating the qualitative difference between sensitive networks ($G_2$, $G_7$, and $G_{19}$ from Fig. 7) and nonsensitive networks ($G_1$ from Fig. 7). The red cross symbol indicates the value of $\mu$ at $\delta =0$ corresponding to the case of no perturbation, which is the same for all four networks. The inset shows a zoom-in plot of the marked rectangular region surrounding the red cross. The perturbation matrix $\mathrm{\Delta }A$ was chosen randomly following the same procedure used in Fig. 7. To facilitate visualization, in this figure we allow negative $\delta$, which corresponds to considering a perturbation term of the form $|\delta |(-\mathrm{\Delta }A)$. (b) Log-log plot of the change in convergence rate $|\mathrm{\Delta }\mu |$ versus $\delta$, which confirms the scaling $|\mathrm{\Delta }\mu |\sim {\delta }^{1/\beta }$ for small $\delta$.