|Home | About | Journals | Submit | Contact Us | Français|
We propose an analytical technique to study large fluctuations and switching from internal noise in complex networks. Using order-disorder kinetics as a generic example, we construct and analyze the most probable, or optimal path of fluctuations from one ordered state to another in real and synthetic networks. The method allows us to compute the distribution of large fluctuations and the time scale associated with switching between ordered states for networks consistent with mean-field assumptions. In general, we quantify how network heterogeneity influences the scaling patterns and probabilities of fluctuations. For instance, we find that the probability of a large fluctuation near an order-disorder transition decreases exponentially with the participation ratio of a network’s principle eigenvector – measuring how many nodes effectively contribute to an ordered state. Finally, the proposed theory is used to answer how and where a network should be targeted in order to optimize the time needed to observe a switch.
Network science is highly interdisciplinary, and when combined with other fields such as statistical physics and nonlinear dynamics, provides a useful framework to address fundamental questions regarding complex systems1–3. Network approaches have provided quantifiable results in diverse applications in nearly every field of science and engineering, from biological networks4 to climate networks5, information networks6, infrastructure networks7, and social networks8. Consequently, much progress has been made in understanding the role of topology in many collective processes in complex systems, including in adaptive and co-evolving networks, where the topological dynamics is itself a function of a network’s state9. Popular examples where interaction structure is understood to strongly influence behavior are the spread of infectious diseases10, the dynamics of neural systems11, the synchronization of coupled oscillators12, the patterns of voters13, and the collective motion of networked swarms14.
However, many theoretical results in network dynamics rely on deterministic, and “mean field” limits of some simple model2. Though useful, such approaches typically ignore noise and dynamical fluctuations that are inherent in virtually all of the aforementioned examples. Although inherent noise may be considered small in large networks, the existence and observation of large fluctuations can result in drastic change in a network’s dynamics15. Therefore, some recent efforts have been made to study the relationship between network dynamics and noise16–18. It has been demonstrated that the interplay between complex topology and noise can alter well-known scaling laws and patterns for fluctuations19, as well as provide new control mechanisms that take advantage of noise-induced phenomena including switching and extinction15, 20.
An important class of statistical physics models for capturing many processes, and where noise is relevant, are spin systems, in which nodes in a network take on discrete states1, 2. Typically, the probability (or probability of changing in time) of any configuration of states depends on the configuration’s energy21. Such an approach has been useful for understanding opinion formation and dynamics22, 23, rumor spreading24, as well for machine learning on networks (e. g., Boltzmann machines)25 and network inference26. In spin dynamics, often two limits are considered. The first entails zero temperature, where the energy of a network is minimized and generally does not fluctuate between configurations22, 27. The second involves the average state of a network at finite temperatures and noise, where an order-disorder transition is observed and analyzed in the limit where the network size tends to infinity23. However, all real networks are finite systems. Thus, fluctuations between distinct metastable configurations arise – effectively changing the collective order due to noise. Example systems are social networks, where fluctuations are known to be important28, 29, and hence switching from one majority opinion to another is possible. Yet, many open questions remain about how switching occurs in complex networks as a result of random fluctuations.
On the other hand, the role of noise is reasonably well understood in simple well-mixed and spatially homogenous contexts30–33. It has been demonstrated in many works that noise and collective dynamics can couple in such a way as to induce a large fluctuation – effectively driving a system to switch from one collective behavior to another. If the fluctuation is a rare event, then the process is captured by a most probable, or optimal path (OP) – where all others are exponentially less likely to occur34, 35. In such cases, the OP is describable in an analytical-mechanics formalism familiar from classical physics30, which has been used to elegantly describe a variety of rare phenomena including: fixation in evolutionary games36, extinction of disease in homogenous populations37, switching in self-regulating genes38, large velocity fluctuations in propagating fronts39, viral clearance40, irreversible fluctuations in electronic circuits41, and switching in quantum mechanical oscillators42. Our strategy is to find the OP in general network configurations given a general spin (opinion) dynamics, and use it to understand the dynamical switching pathway between ordered (majority) states, the average time needed to see a switch, and the distribution of large fluctuations in both real and synthetic networks.
This report moves well beyond both the deterministic and small fluctuation limits of order-disorder dynamics in complex networks by addressing how large fluctuations occur. Our approach enables us to construct and analyze the approximate OP through a given network – reducing a very high dimensional stochastic and rare process to a single trajectory. Consequently, we are able to compute several topologically dependent quantities that have so far eluded analysis in network science: the distribution, shape, and time scale of large fluctuations. Beyond computation, we show that just above an order-disorder transition, the probability of a large fluctuation decreases exponentially with the participation ratio of a network43, and hence is exponentially sensitive to topological heterogeneity– e.g., the ratio of second to fourth moments of a network’s degree distribution. Moreover, we find two quantitatively distinct scaling patterns for the distribution of large fluctuations, in which low eigenvector-centrality fluctuations predominate at high order (large majorities), and high centrality at low order (small majorities). Finally, we demonstrate with several examples on a Facebook network44 how the formalism is useful for designing controls that optimally leverage noise in order to minimize the time scale for switching – answering where and at what rate a network should be targeted in order to induce a switch.
We consider a system of N nodes interacting through a network. The network is represented by a real-valued matrix, A, where A ij gives the influence strength of node i on node j. At any instant, each node i is in one of two possible opinion (spin) states, characterized by 16, 21. Nodes can change state by interacting with their nearest neighbors in the network, such that s i evolves stochastically in time with a probability that depends on s i and s j of neighbors, for A ij≠0. We study a simple interaction rule motivated from statistical physics, where each node has a tendency to align its opinion with neighboring opinions such that the energy, Ei ≡ − si∑jAijsj, is minimized2, 21, 45. We choose the kinetics to be a continuous-time Glauber dynamics– a Markov process with a transition rate for each node:
where β is an inverse temperature that measures the ratio of energy to thermal noise, f i is a local spontaneous flipping rate16, and α is a rate constant that determines the units of time; we take α=1 without loss of generality. Qualitatively, minimizing energy pushes the network toward complete order (unanimous majority or consensus), while thermal noise and spontaneous flipping inject fluctuations that tend to break up order. Because of this interplay, a majority order emerges as long as β is above some critical value (see supplementary information section two, SI. 2, for derivation)23, 46. Below, our results are to be compared with the stochastic process defined by Eq. (1).
In general, the network dynamics is governed by a master equation for the N-node probability distribution, ρ(s, t), which is high-dimensional and difficult to analyze in its entirety. Therefore, we seek a simplified description of the dynamics by first considering the opinion density at each node in a large ensemble. The ensemble consists of C identical networks with the same A, but independent realizations of the dynamics given by Eq. (1). The opinion density in the ensemble is , where s i,c is the state of node i in realization c of the stochastic process. Ultimately, we are interested in the limit C→∞, or continuous density. Our goal is to find an approximate master equation for the network ensemble distribution, P(m, t), that is a function of the densities alone, and extract a particular solution relevant to large fluctuations. As we will see, such a solution will correspond to the OP. To find it, we must consider the transition rates for m i and make a mean-field approximation:
The mean-field approximation replaces s i,c by its ensemble average m i – effectively neglecting correlations between neighbors in the network. The result is a master equation that describes a simplified stochastic process in terms of the opinion density at each node,
When β is above threshold, P(m, t) is peaked around one of two ordered equilibrium states, m(t)≈±m*. Dynamically, a finite network fluctuates around equilibrium, m*, after an initial transient, for a long time until a large fluctuation occurs, which carries the network to the opposite ordered state, −m *, as shown in Fig. 1(a). Such order switches are rare events in large networks (N ≫ 1), and we expect them to be encoded in the tails of P(m, t). In particular, if m corresponds to a large deviation from m*, we expect an exponential reduction in probability, as demonstrated in Fig. 1(b). Therefore, it is convenient to constrain our search for solutions of the master equation to the exponentially-distributed tail that is relevant for rare events, since Eq. (3) contains too much information to be useful in practice.
These observations suggest extracting a solution of Eq. (3) with exponential, or Wentzel–Kramers–Brillouin (WKB) form, P(m, t)=ae −CS(m,t). The WKB solution for the ensemble distribution can be viewed as a product of independent and identical distributions for each realization in the ensemble. Thus, we can approximate the probability distribution for states in a single realization, ρ(s, t), by
We can find the continuous-density solution for P(m, t) in the large-ensemble limit by substituting the WKB ansatz into Eq. (3), Taylor expanding Eq. (3) in powers of the small parameter 1/C, and neglecting terms of or smaller (see SI. 1 for expansion). As is customary, this approximation converts the master equation into a familiar Hamilton-Jacobi equation (HJE) from analytical mechanics30, 33:
where S and H are called the Action and Hamiltonian, respectively33. Once S is found from Eq. (5), the distribution of large fluctuations is determined. Following analytical mechanics, we define a momentum, pi = ∂S/∂mi, and conveniently express the Hamiltonian,
A crucial result of the WKB approximation, which makes the approach worthwhile, is that solutions of the HJE extremize S, when expressed as the integral30
where m(t) and p(t) are determined from Hamilton’s equations of motion below, Eqs (8–9). Since S is minimized, the probability of a stochastic path associated with S is maximized. In summary, when the transition between ordered states in a network is exponentially rare, there is a least-action path in (m, p) phase-space that connects the two, which is a local maximum in probability, and therefore corresponds to the OP through a network. The Action along the OP gives ρ(m, t) from Eqs (4) and (7).
Finally, just as in analytical mechanics, a convenient approach for computing the OP is to solve Hamilton’s equations of motion for the system, and :
It is important to notice that if one takes p ≡ 0 in Eqs (8–9) the “quenched mean field” equations are derived for a kinetic Ising model as a special case which ignores fluctuations46, 47. Similar findings have been shown recently for epidemic dynamics48. Therefore, the OP formalism naturally generalizes deterministic approaches for network dynamics to include large fluctuations. In practice, all that is needed to find the OP are appropriate boundary conditions for Eqs (8–9). These are derivable directly from the distribution as explained in Sec. 3.2.
By considering histograms of m(t) from time-series data, e.g., Fig. 1(a), we can see that the expected exponential distribution of large fluctuations appears. Moreover, simple inspection gives us the boundary conditions for solving Eqs (8–9). Since the distribution takes a maximum value at the equilibrium m* (satisfying ), the initial boundary condition for a large fluctuation is m(t=0)=m* and . Because we are interested in large fluctuations that lead to a switch, the final boundary condition is similarly m(t→∞)=−m * and p(t→∞)=0 30, 33, 49. Therefore, by solving Eqs (8–9) subject to zero-momentum boundary conditions, we determine the OP: a single trajectory that gives the probability of a large fluctuation to m within logarithmic accuracy, ln ρ(m) ≈ − S(m). We note that since the distribution is nearly constant in time, the network Action is time-independent ∂S/∂t = 0 = H ∀t. Therefore, S(m) is equal to the line integral of the momentum along the OP, from Eq. (7).
Optimal paths can be computed numerically from Eqs (8 and 9) with zero-momentum boundary conditions using quasi-Newton methods49 and dimension-reduction techniques based on a spectral decomposition of A (see SI. 6–7 for details on numerical approaches). In general, we find two distinct segments of the OP: an activation segment with p≤0, requiring noise to carry the network from m=m* to m=0 (shown in Fig. 3), and a deterministic segment with p=0 that leads from m=0 to m=−m *. Since the deterministic segment has p=0, the probability distribution is predicted to be flat from m=0 to m=−m*, as demonstrated in Fig. 1(b). In addition to the distribution, the shapes of large fluctuations are predicted and can be compared to simulations. For example, Fig. 2 shows OP projections and prehistory trajectory-density plots from many stochastic realizations of Eq. (1) for two networks. We see that the OP for each network lies near the maximum of the corresponding heat maps. Importantly, since the OP predicts the most likely sequence of changes in opinion density, m-projections at various points along a path (such as in Fig. 2) give us insight into how and where noise acts in a network during a large fluctuation.
Beyond finding a prescription for computing the OP numerically, we are interested in understanding how the path structure is related to topological properties of a network. The OP dependencies can be found analytically in certain limiting cases of interest, such as near threshold. A general method used throughout Sec. 3.2 is to expand Eqs (8 and 9) around the OP boundaries – giving the local functional forms of m, p, and S(m). We show in SI. 4 that when f=0 and the distance to threshold is small, δ ≡ βλ − 1 ≳ 0 (where λ is the largest eigenvalue of A), the OP depends on the principle right and left eigenvectors of A, η and ζ, respectively:
h is a unit-length parameter along the activation segment, . In general, the principle right eigenvector satisfies η=A η/λ, and similarly ζ=ζ A/λ with ζ·η=1; η i is called the right eigenvector centrality of node i, and is an approximate measure of the importance of node i in the network3. Near threshold m and p are parallel to the eigenvectors of A. Therefore, if η and ζ contain relatively few nodes that have significantly large eigenvector centrality compared to most others, we expect the opinion density and fluctuations to be similarly large at such nodes. Further examining Eqs (10 and 11), we see that the momentum is a cubic function of h to lowest order in δ, implying that ln ρ(m) has a curvature with respect to opinion density that is twice as large at m* compared to m=0, and therefore the probability changes more quickly near m* than near m=0. This can be contrasted with large fluctuations that cause extinction in epidemics, for which the momentum is linear15.
In general, since each node’s contribution to the Action is equal to the line integral of the momentum, we expect the contribution to increase with increasing eigenvector centrality. This pattern can be seen in Fig. 3, where the area under the darker curves (corresponding to higher centrality) contains the area under the lighter curves. Near threshold, we can calculate the line integral explicitly by substituting Eqs (10 and 11) into Eq. (7): or
The expression is interesting, since for a symmetric network we find that S is a function of the fourth moment of the eigenvector-centrality distribution, or inverse participation ratio43. To understand what this implies we first consider the case of a homogeneous complete graph, where In this case the probability of a large fluctuation to zero consensus, and therefore the probability of switching, is to logarithmic accuracy (for Nδ2 ≫ 1). On the other hand, for a random network without degree correlations the standard annealed-network approximation gives given a degree k i for node i and a network average of k 2, k 2. Hence, . If the the network is composed of nodes with only degree k, the switching probability is predicted to equal the complete graph’s value. However for a heterogeneous network, e.g., when the distribution of nodes with degree k scales like k −γ, the additional topological factor, 〈k2〉2/〈k4〉, can be significantly smaller than 1 depending on the cutoff of the distribution. Generally then, given the same distance to threshold and system size the switching probability is predicted to increase exponentially with increasing topological heterogeneity15. This will be particularly the case if γ≤5, which can be contrasted to other kinetic systems such as epidemics48.
Near threshold a network is highly disordered, and we may wonder how large fluctuations behave as the order increases with β far above threshold. In fact, the scaling patterns for the very largest fluctuations to small m≈0 maintain a similar form. In SI. 5 we show that m i and p i remain proportional to centrality in the tail of the distribution by expanding Eqs (8 and 9) around the origin (assuming f i=f and A is symmetric). Hence, the relative probabilities for observing small ordering are simple functions of a network’s average-squared opinion density,
The generic scaling with m 2 occurs because p i tends to the same slope with respect to m i along the OP for all nodes when m≈0; the scaling is demonstrated in Fig. 3 with a black-dashed line. Notably Eq. (13) implies that the tail of ρ(m) is an exponential with a rate that increases linearly with the network size and distance to threshold, βλ−1−2f, but decreases with f. The latter effect entails an additional broadening of the distribution in the presence of spontaneous flipping that is independent of network topology.
In contrast, the scaling at high order is significantly different, with fluctuations that are very sensitive to node centrality. In particular, when a network is near consensus, m ≈ m⁎ ≲ 1, lower-centrality nodes have fluctuations that are exponentially larger than higher-centrality nodes along the OP, and are distributed according to a simple Gaussian:
Intuitively, when a network is near consensus, it is very improbable for high centrality nodes to change state. This pattern can be seen in Fig. 3 where the slope of p i with respect to m i along the OP is much steeper for higher centrality nodes when m≈m*– corresponding to a more rapid decrease in probability from the maximum value. The quantitative scaling for high order is depicted in Fig. 3 with a blue-dashed line, and reflects the pattern that the standard deviation of each node’s distribution is exponentially decreasing with its centrality near consensus. Equation (14) results from an expansion of Eqs (8 and 9) near consensus, given in SI. 5, and assumes f≈0 and A≈λ ηη T.
The final observable we consider in this report is the average switching time, T, since it quantifies the expected time scale over which the very largest fluctuations occur. Generally T takes the form:
from the assumption that switching between ordered states has a rate, or inverse time, proportional to the probability, ρ(0) 33, 34. For sufficiently large S, the exponential contribution dominates, and therefore lnT~S(0), as demonstrated in Fig. 4(a) for several networks. Given the dependence on ρ(0), our analysis indicates that switching times have exponential sensitivity to the topological properties which determine the Action, such as network heterogeneity.
Since we are able to predict large fluctuations with the OP theory, we may be interested in adding controls to the network dynamics in order to optimize an observable– e.g., minimize T. In this report, our approach will involve minimizing S(0) (a deterministic quantity derived from theory). However, it is important to note that the formalism ultimately entails using internal fluctuations in the network to do work, which would not be possible in the absence of noise15. In this section we study the example of targeting a subset, F of nodes with spontaneous flipping, for which f i=f≠0 given iF, while all other nodes have rate equal to zero. We consider several cases using the Facebook network44 as an illustration, and answer the following questions: what nodes should be targeted in order to minimize T (given a fixed number for targeting, |F|, and a fixed flipping rate, f), and whether targeting a larger subset of nodes at a lower rate, or a smaller set with a higher rate tends to minimize T.
Following the analysis above, we differentiate among nodes by eigenvector centrality, η i for node i. In order to first compare equal numbers of nodes for targeting, we have arranged nodes in a list according to increasing η i, and binned the list into roughly equally sized bins. In Fig. 4(b) ln T is shown as a function of the average η i in F, ηF, for several f. The size of F is fixed at 32 nodes. The F with highest ηF in Fig. 4(b) represents nodes near the maximum value of ηF in the network: i.e., the 32 nodes with highest centrality. The second highest ηF represents nodes with the next highest centrality, and so on (see SI. 7–8 for more details). Since the times decrease with ηF, we can see that it is optimal to target nodes with higher centrality. The finding makes intuitive sense, since high η i implies high importance for a given node3, and therefore one might expect increased effect from flipping high-centrality nodes. Nevertheless, the scale of difference is important. Because of the exponential form of T, targeting nodes with the highest centrality can reduce lnT by 25% for the parameters shown, even though <1% of the network is controlled. Predictions are in good agreement with simulations, Fig. 4(b).
On the other hand, with a different control scheme targeting higher η i alone may be sub-optimal. Another approach is to start with the 32 highest-centrality nodes (i.e., the control with the smallest lnT in Fig. 4(b), and increase/decrease the size of F by adding/removing nodes with lower η i. An example is shown in Fig. 4(c) (blue circles) where |F| nodes with the highest η i are targeted. In order to keep the total rate of flipping constant, f|F| – a proxy for the amount of work done by the controller, f must vary accordingly. We can see that for the Facebook network it is more optimal to target a larger set of nodes at a slower rate, than fewer nodes at a higher rate, even though a larger set implies decreasing ηF. Similar to Fig. 4(b), a reduction by a factor of 25% in lnT is predicted and observed, despite targeting at most <2% of the network. Finally, it is interesting that minimizing S(0), with the controls considered, is strongly correlated with reducing the amount of order in the network. To demonstrate, we redo the second control by picking an f so that m *2 is constant as we vary F. The result is a slowly varying T, illustrated in Fig. 4(c) (green diamonds).
There is much interest in understanding the relationship between dynamics in complex systems and their underlying topology. However, most theoretical results pertain to deterministic limits, where noise is ignored. Therefore, analytical and computational tools are needed to understand how noise and dynamics interact, especially when the interplay causes a large qualitative change in a network’s collective dynamics. In this report, we have developed an approach based on WKB techniques applied to finite networks, which allowed us to analyze large fluctuations in order-disorder transitions driven by internal noise. The work went well beyond steady state and threshold analysis – concerning itself with a global dynamical object, the optimal path of large fluctuations through a network. By computing the optimal path we were able to predict the probabilities and time scales associated with large fluctuations, including the largest fluctuations that entailed a noise-induced switch between deterministically stable states. Simulations in both real and synthetic networks showed good agreement with theory.
The optimal path approach applied to networks has several advantages, which were detailed in this report. First, it provides a way to predict large fluctuations that naturally reproduces popular mean-field results as a special zero-momentum limit. Because of this structure, we expect more accurate network-approximation techniques, such as those that include dynamical correlations, to be generalized in a similar way – and for many other dynamical processes in networks. For instance, WKB techniques can be used to analyze large fluctuations in nonlinear systems with time delays50, colored noise51, and memory effects52. Second, it allows one to analytically quantify the scaling patterns of large fluctuations on topology, including the exponential sensitivity of fluctuation-probabilities to topological heterogeneity and the multi-step structure of large fluctuations from highly ordered states through heterogeneous networks. Third, optimal control of noisy network dynamics is reduced to deterministic control of a mechanical analog. By minimizing the network Action derived from the optimal path theory, one can construct the optimal way to leverage internal noise so as to maximize the probability of network switching. This was demonstrated in a Facebook network where minimum-Action controls correctly predicted large exponential reductions in the average switching time, by targeting an optimal subset of the network that represented less than two percent of the total. Fourth, we expect that our approach will be useful for current avenues of research in network science and many new applications, such as network inference from data in the presence of large fluctuations.
J.H. is a National Research Council postdoctoral fellows. I.B.S. was supported by the U.S. Naval Research Laboratory funding (N0001414WX00023) and office of Naval Research (N0001416WX00657) and (N0001416WX01643). We are very grateful to L.B. Shaw and L. Mier-y-Teran-Romero for useful discussions.
Both J.H. and I.B.S. contributed to the methods and analysis, and wrote the manuscript.
The authors declare that they have no competing interests.
Electronic supplementary material
Supplementary information accompanies this paper at doi:10.1038/s41598-017-08828-8
Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.