key: cord-0010922-03iukrfw
authors: Hindes, Jason; Assaf, Michael
title: Degree Dispersion Increases the Rate of Rare Events in Population Networks
date: 2019-08-09
journal: nan
DOI: 10.1103/physrevlett.123.068301
sha: ba98dd5bc8a33607840b4e370faf1486574f3b41
doc_id: 10922
cord_uid: 03iukrfw

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] [3] [4] [5] [6] , switching in gene regulatory networks [7] [8] [9] [10] , the arrival of biomolecules at small cellular receptors [11] , and power-grid destabilization [12] [13] [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] [17] [18] [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] [21] [22] , often the resulting Hamilton equations are difficult to solve, as they require computing unstable trajectories in high-dimensional phase spaces [23] [24] [25] [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 (S) and infecteds (I) [27] . A susceptible can get infected upon encountering an infected individual, S þ I → I þ I, while an infected can recover and become susceptible again, I → S. We first consider networks with only two degree classes, and then generalize to arbitrary degree distributions. We assume a network of N ≫ 1 nodes, with N=2 nodes of degree k 1 ≡ k 0 ð1 − ϵÞ and N=2 nodes of degree k 2 ≡ k 0 ð1 þ ϵÞ. 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 n i the number of degree-k i (i ¼ 1, 2) infected nodes, and by x i ¼ n i =ðN=2Þ the densities of degree-k i 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 Φðn 1 ; n 2 Þ≡ Φðx 1 ; x 2 Þ ¼ ðk 1 x 1 þ k 2 x 2 Þ=ðk 1 þ k 2 Þ. Thus, the average infection rate (per individual) of a susceptible node of degree k i is λk i ð1 − x i ÞΦðx 1 ; x 2 Þ, while the recovery rate is simply x i .

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] [19] [20] [21] [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 hki, allowing for a rigorous perturbative treatment. For bimodal networks hki ≡ k 0 , while σ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi

The deterministic rate equations, describing the mean density of infected nodes with degrees k 1 and k 2 , read

The critical value of λ, below which there is no longlived endemic state, satisfies on random networks λ c ≡

Þ=k 0 (given the ANA) [28] . Thus, we write λ ¼ Λλ c , where Λ ≥ 1, and Λ − 1 measures the distance to bifurcation or threshold. Rate equations (2) admit two positive fixed points. For ϵ ≪ 1, these become ½x 1 ; x 2 ¼ ½x 0 ð1 − ϵ=ΛÞ;

x 0 ð1 þ ϵ=ΛÞ, which is stable, and ½x 1 ;

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 n 1 ¼ n 2 ¼ 0 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] [31] [32] .

Accounting for demographic noise, the master equation for P n 1 ;n 2 ðtÞ: the probability to find at time t, n 1 and n 2 infected nodes on degrees k 1 and k 2 , respectively, satisfies

where λ ¼ Λð1 − ϵ 2 Þ=k 0 , and E j n fðnÞ ¼ fðn þ jÞ 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 N is large, and the mean time to extinction (MTE), T, 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 1=T, while simultaneously the extinction probability grows and reaches the value of 1 at infinite time [1, 4] . We now plug the ansatz P n 1 ;n 2 ≃ π n 1 ;n 2 e −t=T into master equation (3), where π n 1 ;n 2 is the QSD, and employ the Wentzel-Kramers-Brillouin (WKB) approximation for the QSD, π n 1 ;n 2 ≡ πðx 1 ; x 2 Þ ∼ e −NSðx 1 ;x 2 Þ , where Sðx 1 ; x 2 Þ is the action function [1] . In the leading order in N ≫ 1 we arrive at a stationary Hamilton-Jacobi equation

Once SðxÞ is known by solving Hamilton's equations, so is the MTE, which is proportional to e NSð0;0Þ [4, 21, 30] .

For convenience, let us define new variables

This transformation is canonical since the determinant of the Jacobian ∂ðQ; PÞ=∂ðx; pÞ ¼ 1, where Q ¼ ðu; wÞ, P ¼ ðp u ; p w Þ, x ¼ ðx 1 ; x 2 Þ, and p ¼ ðp 1 ; p 2 Þ. Using the new variables, the path to extinction connects between the fixed points ½w à ; u à ; 0; 0 and ½0; 0; p à w ; p à u , where

Since the above transformation is canonical, the action along the path to extinction is given by [1] Sð0Þ

Transforming to the new variables in Hamiltonian (4), and assuming u and p u scale as OðϵÞ, we find the trajectories p w ðwÞ and p u ðuÞ up to Oðϵ 2 Þ [33] . The trajectories are then substituted into Eq. (6), which yields PHYSICAL REVIEW LETTERS 123, 068301 (2019)

where S 0 ¼ 1=Λ þ ln Λ − 1 is the action for a degreehomogeneous network (ϵ ¼ 0), and f E ðΛÞ > 0. 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 ϵ ≪ 1 our analytical results (7) agree well with numerical solutions of the Hamilton equations, obtained using the iterative action minimization method [26, 33] . 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 gðkÞ the node degree distribution. That is, if N k are the number of nodes of degree k such that P k N k ¼ N, we have gðkÞ ¼ N k =N. We assume that gðkÞ is a symmetric distribution about the mean k 0 ≡ hki, such that gðk 0 þ iÞ ¼ gðk 0 − iÞ for i ¼ 1; 2; 3; …. Let us also assume our distribution has a bounded support such that k min ¼ k 0 − Δ and k max ¼ k 0 þ Δ, where gðk < k min Þ ¼ gðk > k max Þ ¼ 0. We again denote by n k the number of infected individuals on degree-k nodes, and by x k ¼ ½1=gðkÞn k =N ¼ n k =N k the fraction of such infected individuals. Writing down the master equation for P fn k g -the joint probability to find ðn k min ; …; n k max Þ infected nodes of degree k, and using the above WKB formalism, PðxÞ ∼ e −NSðxÞ , where x ¼ ðx k min ; …; x k max Þ, we arrive at a Hamiltonian equivalent to [21] . Denoting gðkÞp k ¼ ∂S=∂x x , the action can be shown to satisfy [33] Sð0Þ ¼

where we have used the symmetry of gðkÞ about its mean k 0 . Now, since each pair of nodes k 0 AE j for j ∈ ½1; Δ can be viewed as a bimodal network, using Eqs. (6) and (7), the action for such a bimodal network with degrees k 0 − j and k 0 þ j, satisfies ð1=2Þ

Moreover, the node of rank k 0 can be viewed as a bimodal network with ϵ j ¼ 0, such that R p k 0 dx k 0 ¼ S 0 . Therefore, using the fact that P k gðkÞ ¼ 1 and that the variance of gðkÞ satisfies σ 2 ¼ P k ðk − k 0 Þ 2 gðkÞ, the action [Eq. (8) ] and MTE become T ∼ e NSð0Þ ; Sð0Þ

Equation (9) is the first of the main results in this Letter. Namely, for any network, if the CV is small, σ=hki ≪ 1, 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 σ=hki ≫ N −1=2 , 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, f E ðΛÞ, 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 σ 2 =hki 2 and Λ.

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 i flips spin at a rate proportional to 1=½1 þ exp fλΔE i g, where ΔE i is the change in the local pairwise ferromagnetic energy for node i to flip spin, and λ is an inverse temperature [33] . Here, the densities x k are the magnetization of nodes with degree k: the fraction of degree-k nodes with spin (þ) minus those with spin (−). The master equation and Hamiltonian for x can be derived in precisely the same way as the SIS model [37] . The Hamiltonian reads

wherex ¼ P k kgðkÞx k =hki is the degree-weighted mean magnetization, and gðkÞp k ¼ ∂S=∂x k are the momenta.

In contrast to the SIS model, the spin model exhibits three fixed points: x ¼ x à and x ¼ −x à which are stable, and x ¼ 0 which is unstable. The stable fixed points emerge at a pitchfork bifurcation when λ ¼ λ c ≡ hki=hk 2 i. As before, we may denote λ ¼ Λλ c , where Λ ¼ 1 is the bifurcation threshold. In the spin model, demographic noise causes switching between x à and −x à [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 x à coincides up to time reversal, with the fluctuational path from x à to 0 [1] . Once at the unstable point 0, the network can switch to −x à following its deterministic dynamics.

In order to find the switching path, we again use Hamilton's equations gðkÞ_ x k ¼ ∂H=∂p k . The relevant trajectories p k ðxÞ can be found by equating −_ x k j p¼0 ¼ _ x k ðpÞ, where the former represents (minus) the deterministic trajectory. By doing so, the switching path satisfies [33] p k ðxÞ ¼ ð1=2Þ ln ½ð1 þ x k Þ=ð1 − x k Þ − λkx;

and hence the action for switching, Sð0Þ ¼ P k gðkÞ R 0

x à k p k dx k , becomes

Following the same general approach as for the SIS model, we write k ¼ k 0 ð1 þ ϵÞ where ϵ ≡ ðk − k 0 Þ=k 0 . For degree distributions with a small CV, σ=k 0 ≪ 1, we have λ ≈ Λ½1 − hϵ 2 i=k 0 and hϵ 2 i ¼ σ 2 =k 2 0 , as before. In order to evaluate Eq. (11) in the limit of hjϵji ≪ 1, we use the small-hjϵji expansion of x à k andx à , see [33] , and keep terms up to order hϵ 2 i. This procedure yields the action and mean switching time (MST)

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.

To check the universality of our results, in Fig. 4 we plot the correction ½Sð0Þ − S 0 =fðΛÞ 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 fðΛÞ, 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] .

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 N −1=2 , 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] . 

Stochastic Population Dynamics in Ecology and Conservation

Modeling Infectious Diseases in Humans and Animals

for more details on the model and additional results from the analysis and simulation

Dynamical Processes on Complex Networks

A Kinetic View of Statistical Physics

In the Supplemental Material [33] we demonstrate that breaking detailed balance still produces a reduced action which is proportional to the coefficient of variation squared