Large order fluctuations, switching, and control in complex networks

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.


Master equation expansion
We want to find the leading contribution to Eq.(3)(main text) when C 1, where C is the number of stochastic realizations of the dynamical process defined by Eq.(1)(main text). Taylor expanding the probability and rates we find and where H(m, p) is given by Eq.(6)(main text). Dividing by C and ae −CS (m,t) , and neglecting the O(1/C) sum, we find Eq.(5)(main text).

Equilibria and linear spectra
In general Eqs.(8-9)(main text) have three equilibria, m = 0 and ± m * , that satisfyṁ =ṗ = 0 with p = 0. In the special case where the spontaneous flipping rate is homogeneous, f i = f ∀i, we find a set of fixed-point conditions: Often, it is convenient to approximate A by the largest term in its eigenvalue decomposition, i.e., a large spectral-gap assumption. When A is symmetric, A = A T ≈ ληη T , where η i is the eigenvector centrality of node i and λ is the largest eigenvalue. Hence, a single equation determines the ordered equilibria (states) in terms of the order-parameter, m * = i η i m * i / j η j : where i η 2 i = 1. The linear stability spectra of the equilibria are found by substituting m = m * + and p = µ into Eqs.(8-9)(main text) and solving to O( ) and O(µ). Because the resulting equations are linear, the dynamics are exponential: (t) = e σt and µ(t) = µe σt . The linearized dynamics gives three equations for the exponent σ and the relative size (shape) of the modes i and µ i : where M is an arbitrary constant.
Ordered states emerge at a threshold where the equilibrium (m = 0, p = 0) changes stability, σ(0) = 0: When βλ > 1 + 2f , ordered sates have σ(m * ) > 0 and σ(0) < 0 for the modes Eqs.(8-9). However, the meanfield assumptions implicit in our approach can be quantitatively inaccurate for the threshold depending on the network. However, the WKB approach can be augmented to include pairwise correlations, for example, which generally improves accuracy. We mention that other solutions are possible with µ i ≡ 0, which have oppositely signed spectra in terms of stability, i.e., where m = 0 is unstable and m = m * is stable. In general, taking p ≡ 0 ∀t in Eqs.(8-9)(main text) gives the so called "quenched mean field" equations corresponding to a dynamic Ising model with random flipping. In this way, the condition p = 0 is what allows a trajectory (i.e, the OP) to exist from a "stable" to "unstable" state in the deterministic mean-field theory.

Network details
The Facebook network used throughout the paper was taken from http://snap.stanford.edu/data/egonets-Facebook.html. It contains 4039 nodes and 88234 edges. The power-law network in Fig.2(a)(main text) was generated from the configuration model with degree(k) distribution, g k = k −2.5 / 300 k =10 k −2.5 , and 600 nodes. In Fig.4(a)(main text) the Erdős-Rényi network had 500 nodes and 15000 edges and the bimodal network was generated from the configuration model with 400 nodes and two degree classes: 40 nodes had degree 50 and 360 nodes had degree 5.

Near threshold OP
As in the main text, we consider the special case where f = 0 near threshold, with δ = βλ − 1 0. Our approach is to find m * and the linear dynamics (Eqs.(7-9)) near m ≈ 0 and m ≈ m * to lowest order in δ. This will give us boundary conditions which can be used to determine the OP to the same order in δ. Once the OP is found, we can explicitly perform the line integral of momentum that gives the Action, Eq.(7)(main text), and hence the probability exponent in the distribution of large fluctuations, Eq.(4)(main text). For this section, it is not assumed that A has a large spectral gap nor is symmetric.
Next, we find the dynamics near m ≈ 0. Analogous to Eqs.(7-9) (without the symmetric A assumption), Solving Eqs.(14-15) we find that i = Eη i and µ i = −δEζ i , which means that or the derivative at the boundary m = 0 and p = 0. Similarly, we seek the derivative at the boundary m = m * and p = 0 to the same order, O(δ). The linearized dynamics are Solving Eqs.(17-18) gives i = Aη i and µ i = 2δAζ i , implying a derivative boundary condition Finally, it is convenient to parameterize m i and p i in terms of a unit-length parameter h, such that h ≡ m i /m * i ∀i. Note: when h = 1 the network is ordered at m * , and when h = 0 the network has no order m = 0. Therefore we can write, is an unknown function that we must determine. The boundary conditions above imply: f (h = 1) = 0, f (h = 0) = 0, df dh (h = 1) = 2, and df dh (h = 0) = −1. If we assume that f (h) is a polynomial, the simplest polynomial that satisfies the four boundary conditions is a cubic function, f (h) = h h − 1 h + 1 . Hence we arrive at Eqs.(10-12)(main text). We mention that the near threshold OP is a convenient initial guess for the Iterative-Action-Minimization-Method described in Sec.6.

Scaling away from threshold
We would like to use Eqs.(7-9) to find basic scalings of the OP away from threshold, where the network's metastable order is high, which will help us understand fluctuations near m ≈ 0 and m ≈ m * -i.e., the largest and smallest fluctuations. We first study the former for which the solution of Eqs.(7-9) is useful: . Therefore, the momentum is linear in m i with constant slope across the network, p i ≈ m i [1 + 2f − βλ]/[1 + 2f ]. By considering the action at m = 0, Using Eq. (19) -showing that fluctuations for low-centrality nodes are exponentially larger than for high-centrality nodes. Moreover, combining with Eq.(7-8) we find σ(m * ) ≈ 1 and which is equivalent to Eq.(14)(main text).

Finding the OP numerically
In general, one would like to find the OP beyond the limiting cases. Of course, no analytic solution is possible except in networks that are effectively low-dimensional. Since the path connects two equilibria via a heteroclinic orbit, in practice it must be constructed numerically, e.g., through shooting, or quasi-newton methods, etc. The method used in this report is of the latter form, namely the Iterative-Action-Minimizing-Method (IAMM) (doi:10.1016/j.phys d.2013.04.001). In this method, OPs are generated from a least-squares algorithm that minimizes the residuals between Eqs.(8-9)(main text) and finite-difference approximations, with fixed-point boundary conditions from Sec.(3.1)(main text) (used to close the differencing). However the dimension for the minimization is 2N d where d is the number of discrete points in the differencing and N is the size of the network, which is prohibitively large for large N (typically we choose 200 ≤ d ≤ 500). Therefore, in practice it is necessary to coarse-grain the network in some way. We describe our approach for this report in Sec.7. We mention that for the special case of f = 0, the OP is reversible, and therefore dm/dt along the first segment is time reversed along the second.

Binning the network
We are interested in reducing the dimension of network defined by the adjacency, matrix, A ∈ L(R N , R N ).
All of the networks considered in this report are symmetric, though the formalism does not require this assumption. Nevertheless, in this section we assume A = A T . Given a sequence of of eigenvalues, {λ j }, and eigenvectors, {η j }, for A, we assume the largest eigenvalue is much greater than all of the others. This is a good approximation for many networks, including those discussed in Sec.3. From the spectral decomposition theorem, we can approximate the adjacency matrix as A ≈ ληη T , where λ = max{λ i }, and η the corresponding eigenvector.
In order to create a mapping from N dimensions to one that is significantly lower, we first notice that the entries of the eigenvector (which we assume is normalized) roughly relate a measure of vertex importance in the graph. Therefore, we reorder the entries of η in increasing order such that v = P η, where P ∈ L(R N , R N ) is a permutation matrix, and v 1 ≤ v 2 ≤ · · · ≤ v N . Notice we have not changed the norm of v, nor have we made any dimension reduction.
Next, we arbitrarily pick a binning of v such that there exists |B| N bins, and associated with each bin b ∈ B we have a distribution, g b , as well as an index set, I b . We define an indicator function on an index such χ b (i) = 1 if i ∈ I b , 0 otherwise. We now define a vector that averages the nodes within a bin b as the following: That is, N g b is the total number of nodes in bin b, and The map in Eq.23 computes the average over all of those nodes in bin b. To compute the entire transformation from R N into R B ,we have where A ∈ L(R N , R B ). Using the same transformation map for m and p, we find the corresponding |B| dimensional vectors, M and P, respectively, for the average opinion density and momentum in bins. By replacing v i with r b , m i with M b , and p i with P b for i ∈ I b in Eqs.(8-9)(main text), we get the (approximate) equations of motion for bin b: assuming A final requirement is needed to ensure that the binned and original system have the same bifurcation point and are similarly normalized: after binning we renormalize r b so that j η 2 j = b r 2 b g b N = 1. In practice, to use the binning procedure one must specify χ b (i). We illustrate with the Facebook network, where we chose |B| = 50. The v i distribution is shown in Fig.1 in blue. Note, the first 3000 nodes have small v i , and therefore we truncate the x-axis for easier viewing. Visually the v i has roughly three relevant parts: each of the three parts with roughly equal numbers of nodes in each bin -with a total of 28, 12, and 10 bins for the three parts, respectively. Such a choice gave the following indicator functions, χ b (i), which we list in