Degree Dispersion Increases the Rate of Rare Events in Population Networks.

There is great interest in predicting rare and extreme events in complex systems, and in particular, understanding the role of network topology in facilitating such events. In this Letter, we show that degree dispersion-the fact that the number of local connections in networks varies broadly-increases the probability of large, rare fluctuations in population networks generically. We perform explicit calculations for two canonical and distinct classes of rare events: network extinction and switching. When the distance to threshold is held constant, and hence stochastic effects are fairly compared among networks, we show that there is a universal, exponential increase in the rate of rare events proportional to the variance of a network's degree distribution over its mean squared.

Systems containing a large, yet finite, population of interacting individuals or dynamical units often experience fluctuations due to the stochastic nature of agent interactions and local dynamics. Most of the time such systems reside in the vicinity of some attractor, undergoing small random excursions around it. Yet, occasionally a rare large fluctuation, on the order of the typical system size, may occur, which can lead to a transition to an absorbing state (a state that, once entered, cannot be left) or to the vicinity of another attractor. As a result, stochasticity can turn deterministically stable attractors into metastable states [1]. Examples of such extreme, rare events, which may be of key practical importance include population extinction [2-6], switching in gene regulatory networks [7-10], the arrival of biomolecules at small cellular receptors [11], and power-grid destabilization [12-14].

Usually, rare events in populations are considered within well-mixed or homogeneous settings, e.g., where individuals interact with an equal number of neighbors. In this case, analytical treatment is possible using standard techniques [6,9,15]. On the other hand, it is known that in topologically heterogeneous networks, e.g., where nodes have a variable degree, the critical behavior can be dramatically affected [16-19]. Unfortunately, predicting rare events in degree-heterogenous networks is notoriously difficult, due to high dimensionality and complex coupling between degrees of freedom. Though some progress has been made by applying semiclassical approximations to master equations governing stochastic dynamics in complex systems [20-22], often the resulting Hamilton equations are difficult to solve, as they require computing unstable trajectories in high-dimensional phase spaces [23-26]. Consequently, analyzing rare events in general networks has been mainly limited to near-bifurcation regimes, where dimensionality is reduced.

In this Letter, we apply a novel perturbation scheme that allows us to predict a universal increase in the rate of rare events by exploiting the extent of network heterogeneity, or degree dispersion. We find that this increase is proportional to the ratio of the variance of a network’s degree distribution to its mean squared, or coefficient of variation (CV) squared, and is otherwise independent of topology. Our approach is shown analytically for two canonical examples of fluctuation-driven rare events: extinction of epidemics in the susceptible-infected-susceptible (SIS) model on networks, and switching (or spontaneous magnetization flipping) in binary spin networks.

Extinction in heterogenous networks: The SIS model.—We begin by considering the SIS model of epidemics, which consists of two types of individuals: susceptibles () and infecteds () [27]. A susceptible can get infected upon encountering an infected individual, , while an infected can recover and become susceptible again, . We first consider networks with only two degree classes, and then generalize to arbitrary degree distributions. We assume a network of nodes, with nodes of degree and nodes of degree . Each node represents a single individual which can be in either state. We assume the infection rate is and the recovery rate is 1.

Denoting by the number of (, 2) infected nodes, and by the densities of infected nodes, the probability for a given node to be connected to an infected node in a random network with this bimodal degree distribution is . Thus, the average infection rate (per individual) of a susceptible node of degree is , while the recovery rate is simply .

In order to make analytical progress, we assume that the average dynamics over an ensemble of uncorrelated random networks can be approximated by the following four (twice the number of degree classes) stochastic reactions, occurring in a well-mixed setting [18-22]: This formulation is equivalent to the so called annealed network approximation (ANA) [28]. However, an analogous argument can be developed for networks with empirical adjacency matrices in the limit of large spectral gaps [29]. In the latter case, the degree is replaced by the eigenvector centrality in all results below.

We are interested in quantifying how broadening a network’s degree distribution affects the rate of extinction of infection by stochastic fluctuations. We focus on the case where the standard deviation of the degree distribution, , is sufficiently smaller than its mean , allowing for a rigorous perturbative treatment. For bimodal networks , while . Therefore, we assume henceforth that , or .

The deterministic rate equations, describing the mean density of infected nodes with degrees and , read The critical value of , below which there is no long-lived endemic state, satisfies on random networks (given the ANA) [28]. Thus, we write , where , and measures the distance to bifurcation or threshold.

Rate equations (2) admit two positive fixed points. For , these become , which is stable, and , which is unstable, where . A transcritical bifurcation occurs as passes the value of 1. While it gives some intuition, the deterministic picture ignores demographic noise emanating from the discreteness of individuals and stochasticity of the reactions. This noise, and the fact that the extinct state is absorbing, make the nontrivial stable fixed point in the language of the rate equations, metastable. Thus, the network ultimately goes extinct via a rare, large fluctuation [4,6,30-32].

Accounting for demographic noise, the master equation for : the probability to find at time , and infected nodes on degrees and , respectively, satisfies where , and is a step operator. Next, we assume that the network settles into a long-lived metastable state prior to extinction. This assumption is justified if is large, and the mean time to extinction (MTE), , is very long. This metastable state, which is described by a quasistationary distribution (QSD) about the stable fixed point, slowly decays in time at a rate which equals , while simultaneously the extinction probability grows and reaches the value of 1 at infinite time [1,4]. We now plug the ansatz into master equation (3), where is the QSD, and employ the Wentzel-Kramers-Brillouin (WKB) approximation for the QSD, , where is the action function [1]. In the leading order in we arrive at a stationary Hamilton-Jacobi equation , with Hamiltonian where are normalized momenta. The Hamilton equations satisfy and . Once is known by solving Hamilton’s equations, so is the MTE, which is proportional to [4,21,30].

For convenience, let us define new variables , , , and . This transformation is canonical since the determinant of the Jacobian , where , , , and . Using the new variables, the path to extinction connects between the fixed points and , where Since the above transformation is canonical, the action along the path to extinction is given by [1] Transforming to the new variables in Hamiltonian (4), and assuming and scale as , we find the trajectories and up to [33]. The trajectories are then substituted into Eq. (6), which yields where is the action for a degree-homogeneous network (), and . As a consequence, we have obtained an exponential increase in the rate of extinction due to network heterogeneity, which only depends on the CV of the network’s degree distribution. In Fig. 1 we demonstrate that in the limit of our analytical results (7) agree well with numerical solutions of the Hamilton equations, obtained using the iterative action minimization method [26,33].

FIG. 1.

Left panel: versus for bimodal networks. Symbols are numerical solutions of the Hamilton equations for , 2, 2.5, 3, 3.5 (top to bottom); lines are the analytical results (7). Right panel: versus . Symbols are numerical solutions for (see left panel). The curve is the second of Eqs. (7).

Given our analysis for bimodal networks, it is straightforward to generalize to arbitrary, symmetric degree distributions, first, and then to skewed distributions. Let us denote by the node degree distribution. That is, if are the number of nodes of degree such that , we have . We assume that is a symmetric distribution about the mean , such that for . Let us also assume our distribution has a bounded support such that and , where . We again denote by the number of infected individuals on nodes, and by the fraction of such infected individuals. Writing down the master equation for —the joint probability to find infected nodes of degree , and using the above WKB formalism, , where , we arrive at a Hamiltonian equivalent to [21]. Denoting , the action can be shown to satisfy [33] where we have used the symmetry of about its mean . Now, since each pair of nodes for can be viewed as a bimodal network, using Eqs. (6) and (7), the action for such a bimodal network with degrees and , satisfies , where . Moreover, the node of rank can be viewed as a bimodal network with , such that . Therefore, using the fact that and that the variance of satisfies , the action [Eq. (8)] and MTE become

Equation (9) is the first of the main results in this Letter. Namely, for any network, if the CV is small, , the logarithm of the MTE decreases linearly with the square of the CV, compared to the degree-homogenous limit. This indicates that for large networks, for which , the extinction rate is exponentially increased when the population resides on a degree-heterogeneous network, compared with the homogenous case—examples include human contact networks, see, e.g., Refs. [34,35]. Furthermore, while the prefactor for the relative increase of the logarithm of the MTE, , is problem specific, it is independent of the network topology, and is computed for any distance to threshold. Figure 2 shows a comparison between Eq. (9) and Monte Carlo simulations for the MTE in several networks, demonstrating the agreement both in terms of and .

FIG. 2.

Left panel: MTE versus the degree dispersion for several networks; for each point, a mean time is computed from 200 stochastic realizations in a fixed network with a given degree distribution. This is repeated for 20 different network realizations with the same degree distribution and the same number of edges. The log of all such averages is then averaged. Error bars are given by the standard deviation of the latter. Results are shown for uniform (green, , , ), Gaussian (red, , , ), and Gamma (magenta, , , ) distributions. Note that each distribution has one tunable parameter for the variance given a fixed . Right panel: MTE versus the threshold parameter . Results are shown for Erdős-Rényi networks (green, , , ) and (magenta, , , ), and Gaussian distributions (red, , , ). The MTEs and error bars were computed through the same procedure as in the left panel.

Our analysis above required that the network degree distribution be symmetric and bounded. However, even for nonbounded asymmetric distributions the MTE is still given by Eq. (9), as long as such distributions are symmetric in the vicinity of their mean and their skewness is small. In fact, one can show that if these conditions are met, the errors contributed from neglected terms, outside of the symmetrical bulk, are negligible [33]. This is demonstrated in Fig. 2 where we show that theoretical expression (9) agrees well with numerics, also in the case of asymmetric gamma distributions. Moreover, in the Supplemental Material we show that our results even hold for power-law networks when the CV is not too large [33].

Switching in heterogenous networks: The Spin model.—Next, we consider a canonical binary spin system, where nodes are either () or (), instead of infected or susceptible, and make stochastic transitions according to a continuous-time Glauber dynamics [35,36]. Namely, if there is no spontaneous transition (analogous to spontaneous recovery in the SIS model), then each node flips spin at a rate proportional to , where is the change in the local pairwise ferromagnetic energy for node to flip spin, and is an inverse temperature [33]. Here, the densities are the magnetization of nodes with degree : the fraction of nodes with spin () minus those with spin (). The master equation and Hamiltonian for can be derived in precisely the same way as the SIS model [37]. The Hamiltonian reads where is the degree-weighted mean magnetization, and are the momenta.

In contrast to the SIS model, the spin model exhibits three fixed points: and which are stable, and which is unstable. The stable fixed points emerge at a pitchfork bifurcation when . As before, we may denote , where is the bifurcation threshold. In the spin model, demographic noise causes switching between and [38]. In order to find the action for switching, we exploit the fact that there is detailed balance in the absence of spontaneous flipping (though this assumption can be relaxed without qualitatively changing our main result [39]). As a consequence, the deterministic trajectory starting from the vicinity of the unstable point 0 and ending at the stable fixed point coincides up to time reversal, with the fluctuational path from to 0 [1]. Once at the unstable point 0, the network can switch to following its deterministic dynamics.

In order to find the switching path, we again use Hamilton’s equations . The relevant trajectories can be found by equating , where the former represents (minus) the deterministic trajectory. By doing so, the switching path satisfies [33] and hence the action for switching, , becomes

Following the same general approach as for the SIS model, we write where . For degree distributions with a small CV, , we have and , as before. In order to evaluate Eq. (11) in the limit of , we use the expansion of and , see [33], and keep terms up to order . This procedure yields the action and mean switching time (MST) where , is the positive solution of , and .

As was the case for extinction, the action for switching is reduced from the homogeneous network limit by a universal correction, which is a product of the network’s CV squared with a model-dependent (though topologically independent) prefactor. As a consequence, the broader the network degree distribution, the more likely switching is to occur between stable magnetization states, given a constant distance to threshold. Figure 3 shows a comparison between Eq. (12) and Monte Carlo simulations for the MST in several networks, analogous to Fig. 2. As with extinction, the results hold for skewed distributions.

FIG. 3.

MST versus (left) the degree dispersion and (right) the threshold parameter. The same networks were used as in Fig. 2; (left): green, ; red, ; magenta,

To check the universality of our results, in Fig. 4 we plot the correction versus the CV, and obtain a collapse across all networks and all , for both models: network simulations and numerical solutions of the Hamilton equations [33]. As our analysis exemplifies, if the rate of rare events (on log scale) is normalized by the correct process-dependent factor , all networks with the same CV collapse onto the same parabola, given a fixed distance to threshold. Moreover, similar plots and results are shown in the Supplemental Material for power-law networks and continuous-noise analogs for both processes [33].

FIG. 4.

Universal correction to the action for extinction and switching versus the CV; ranges from 1.4 to 3.5. Solid markers denote network simulations and follow Figs. 2 and 3. Numerical computations are shown with open markers for extinction (red) and switching (blue) [33]. Dashed line is .

To conclude, we employed a novel perturbation theory that utilizes the extent of heterogeneity in a network, on two prototypical examples of rare events in networks: extinction in the SIS model of epidemics, and spontaneous magnetization switching in a dynamical spin network. We computed the rate of increase of rare events, and showed that it depends solely on the coefficient of variation (CV) of the network’s degree distribution, but is independent of the exact type of network and connectivity matrix. A key insight therein, was to compare different networks with the same distance to threshold, such that deterministic or fluctuation-free stability was held constant, while propensities for noise-induced fluctuations could be isolated. We found that the rate of extinction or switching can be dramatically increased, as long as the CV of the network’s degree distribution exceeds , which is a reasonable assumption for realistic networks. Finally, we have shown that our approach is valid in processes with maintained as well as broken detailed balance, holds across a broad range of network topologies, and generalizes to different noise sources [33]. Thus, we conjecture that our results are applicable to rare events in a wider range of network processes driven by noise, which include local interactions, and where fluctuations drive a network from a metastable state to an unstable state who merge in a single fixed-point bifurcation [33].