Degree Dispersion Increases the Rate of Rare Events in Population Networks

TL;DR

This paper demonstrates that greater degree dispersion in population networks universally amplifies the likelihood of rare, large fluctuations, such as extinction or switching, by an exponential factor related to degree variance.

Contribution

It provides a theoretical framework linking degree dispersion to increased rates of rare events in population networks, supported by explicit calculations for key event types.

Findings

Degree dispersion exponentially increases rare event rates.

The increase is proportional to the variance over the mean squared of degree distribution.

Results are applicable across different types of rare events in networks.

Abstract

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 work, 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.

Equations

$$n_1\stackrel{\lambda k_1(N/2-n_1)\Phi (x_1,x_2)}{\to }{n}_1+1,n_1\stackrel{n_1}{\to }{n}_1-1,$$

$$n_2\stackrel{\lambda k_2(N/2-n_2)\Phi (x_1,x_2)}{\to }{n}_2+1,n_2\stackrel{n_2}{\to }{n}_2-1.$$

$${\dot{x}}_1$$

$${\dot{x}}_2$$

$${\dot{P}}_{n_1,n_2}(t)$$

$$H(x_1,p_1,x_2,p_2)=\frac{\lambda k_0}{2}\Phi (x_1,x_2)[(1-ϵ)(1-x_1)(e^{p_1}-1)$$

$$+(1+ϵ)(1-x_2)(e^{p_2}-1)]+\frac{x_1}{2}(e^{-p_1}-1)+\frac{x_2}{2}(e^{-p_2}-1),$$

$$w^{\ast }$$

$$p_w^{\ast }$$

$$S(0)=\frac{1}{2}\int p_{1}d{x}_1+\frac{1}{2}\int p_{2}d{x}_2=\frac{1}{2}\int p_{w}dw+\frac{1}{2}\int p_{u}du.$$

$$S(0)$$

$$f_E(\Lambda )$$

$$S(0)=\sum _{k=k_0-\Delta }^{k_0+\Delta }g(k)\int p_{k}d{x}_k=g(k_0)\int p_{k_0}d{x}_{k_0}$$

$$+\sum _{j=1}^{\Delta }g(k_0-j)\int p_{k_0-j}d{x}_{k_0-j}+p_{k_0+j}d{x}_{k_0+j},$$

$$T\sim e^{NS(0)},S(0)=S_0-f_E(\Lambda ){\sigma }^2/{⟨k⟩}^2.$$

$$H(\mathbf{x},\mathbf{p})=$$

$$+$$

$$p_k(\mathbf{x})=(1/2)\ln \left[(1+x_k)/(1-x_k)\right]-\lambda k\bar{x},$$

$$S(\mathbf{0})=\frac{\lambda\!\left<k\right>\!{\bar{x}^{*}}{{}^{2}}}{2}-\frac{1}{2}\!\sum_{k}g(k)\!\Bigg{[}\!\ln\!\big{\{}\!1-{x_{k}^{*}}^{2}\big{\}}+x_{k}^{*}\ln\!\bigg{\{}\!\frac{1+x_{k}^{*}}{1-x_{k}^{*}}\!\bigg{\}}\!\!\Bigg{]}\!.$$

$$T\sim e^{NS(0)};S(0)$$

$$f_S(\Lambda )$$

$$\displaystyle\frac{\partial P({\bf n},t)}{\partial t}=\sum_{k}\Bigg{[}\!(n_{k}+1)P({\bf n}+{\bf 1}_{k},t)-n_{k}P({\bf n},t)\;\;+$$

$$\displaystyle\lambda(N_{k}-n_{k}+1)k\bar{x}P({\bf n}-{\bf 1}_{k},t)-\lambda(N_{k}-n_{k})k\bar{x}P({\bf n},t)\!\Bigg{]}.$$

$$H(\mathbf{x},{\partial }_{\mathbf{x}}S)=$$

$$\displaystyle x_{k}\Big{(}\!\exp\big{\{}\!-\partial_{{x}_{k}}S/g_{k}\big{\}}-1\Big{)}\Bigg{]}.$$

$$S(\mathbf{x})=\sum _k\int {\lambda }_{k}d{x}_k=\sum _k{g}_k\int p_{k}d{x}_k.$$

$$w^{\ast }$$

$$p_w^{\ast }$$

$$S(0)=\frac{1}{2}\int p_{1}d{x}_1+\frac{1}{2}\int p_{2}d{x}_2=\frac{1}{2}\int p_{w}dw+\frac{1}{2}\int p_{u}du.$$

$$H(w,p_w,u,p_u)=2w(e^{p_w/2}-1)[\Lambda (1-w)-e^{-p_w/2}]+\mathcal{O}({ϵ}^2).$$

$$p_w(w)=p_w^{(0)}+\left[3(1-w)-\frac{1+2\Lambda }{{\Lambda }^2}+\frac{w(3+w)}{\Lambda (1-w)}\right]{ϵ}^2.$$

$$p_u(u)=\frac{2(\Lambda -1)ϵ}{\Lambda }\left[1+\frac{{\Lambda }^{2}u}{ϵ(\Lambda -1)}\right],$$

$$\frac{1}{2}\int p_{w}dw$$

$$\frac{1}{2}\int p_{u}du$$

$$S(0)=\sum _{k=1}^{\infty }g(k)\int p_{k}d{x}_k.$$

$$S(0)$$

Full text

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 statesDykmanRev . Examples of such extreme, rare events, which may be of key practical importance include population extinction Lande2003 ; Doering ; Assaf2010 ; Meerson2013 ; Assaf2017 , switching in gene regulatory networks Assaf2011 ; MotterPRX2015 ; Biancalani2015 ; Bressloff2017 , the arrival of biomolecules at small cellular receptors Coombs2009 , and power-grid destabilization Nesti ; TimmePRE2017 ; HindesGrids2019 .

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 Bressloff2017 ; Weber2017 ; Assaf2017 . On the other hand, it is known that in topologically heterogeneous networks, e.g., where nodes have variable degree, the critical behavior can be dramatically affected Moore2000 ; Dorogovtsev2002 ; Sood2005 ; Sood2008 . Unfortunately, predicting rare events in degree-heterogenous networks is notoriously hard, due to high dimensionality and complex coupling between degrees of freedom. Though some progress has been made by applying semi-classical approximations to master equations governing stochastic dynamics in complex systems  Assaf2012 ; HindesPRL2016 ; Sabsovich2017 , often, the resulting Hamilton equations are difficult to solve, as they require computing unstable trajectories in high-dimensional phase spaces WEPRE2002 ; Lindley2 ; Schwartz2 ; Nieddu . 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)Keeling1 . A susceptible can get infected upon encountering an infected individual, $S+I\to I+I$, while an infected can recover and become susceptible again, $I\to S$. We first consider networks with only two degree classes, and then generalize to arbitrary degree distributions. We assume a network of $N\gg 1$ nodes, with $N/2$ nodes of degree $k_{1}\!\equiv\!k_{0}(1-\epsilon)$ and $N/2$ nodes of degree $k_{2}\!\equiv\!k_{0}(1+\epsilon)$. Each node represents a single individual which can be in either state. We assume the infection rate is $\lambda$ 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 $\Phi(n_{1},n_{2})\equiv\Phi(x_{1},x_{2})=(k_{1}x_{1}+k_{2}x_{2})/(k_{1}+k_{2})$. Thus, the infection rate (per individual) of a susceptible node of degree $k_{i}$ is $\lambda k_{i}(1-x_{i})\Phi(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 Sood2005 ; Sood2008 ; Assaf2012 ; HindesPRL2016 ; Sabsovich2017 :

[TABLE]

This formulation is equivalent to the so called annealed network approximation (ANA) Pastor . However, an analogous argument can be developed for networks with empirical adjacency matrices in the limit of large spectral gaps HindesPRE2017 . 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, $\sigma$, is sufficiently smaller than its mean $\langle k\rangle$, allowing for a rigorous perturbative treatment. For bimodal networks $\langle k\rangle\!\equiv\!k_{0}$, while $\sigma\!=\!\sqrt{\langle k^{2}\rangle-\langle k\rangle^{2}}\!=\!k_{0}\epsilon$. Therefore, we assume henceforth that $\sigma\ll\langle k\rangle$, or $\epsilon\ll 1$.

The deterministic rate equations, describing the mean density of infected nodes with degrees $k_{1}$ and $k_{2}$, read

[TABLE]

The critical value of $\lambda$, below which there is no long-lived endemic state, satisfies on random networks $\lambda_{c}\!\equiv\!\langle k\rangle/\langle k^{2}\rangle\!=\!1/[k_{0}(1+\epsilon^{2})]\!\simeq\!(1-\epsilon^{2})/k_{0}$ (given the ANA) Pastor . Thus, we write $\lambda\!=\!\Lambda\lambda_{c}$, where $\Lambda\geq 1$, and $\Lambda\!-\!1$ measures the distance to bifurcation, or threshold.

Rate equations (Degree Dispersion Increases the Rate of Rare Events in Population Networks) admit two positive fixed points. For $\epsilon\!\ll\!1$, these become: $[x_{1},x_{2}]\!=\![x_{0}(1-\epsilon/\Lambda),x_{0}(1+\epsilon/\Lambda)]$, which is stable, and $[x_{1},x_{2}]\!=\![0,0]$, which is unstable, where $x_{0}\!=\!(\Lambda-1)/\Lambda$. A transcritical bifurcation occurs as $\Lambda$ 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 non-trivial stable fixed point in the language of the rate equations, metastable. Thus, the network ultimately goes extinct via a rare, large fluctuation KamenevPRE2008 ; Assaf2010 ; Assaf2017 ; Clancy2018 ; Holme2018 .

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

[TABLE]

where $\lambda\!=\!\Lambda(1-\epsilon^{2})/\!k_{0}$, and $E_{n}^{j}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 (see below). This metastable state, which is described by a quasi-stationary 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 DykmanRev ; Assaf2010 . We now plug the ansatz $P_{n_{1},n_{2}}\simeq\pi_{n_{1},n_{2}}e^{-t/T}$ into master equation (3), where $\pi_{n_{1},n_{2}}$ is the QSD, and employ the WKB approximation for the QSD, $\pi_{n_{1},n_{2}}\equiv\pi(x_{1},x_{2})\sim e^{-NS(x_{1},x_{2})}$, where $S(x_{1},x_{2})$ is the action function DykmanRev . In the leading order in $N\gg 1$ we arrive at a Hamilton-Jacobi equation $H(x_{1},x_{2},\partial_{x_{1}}S,\partial_{x_{2}}S)=0$, with Hamiltonian

[TABLE]

where $p_{i}/2\!=\!\partial_{x_{i}}S$ are normalized momenta. The Hamilton equations satisfy $\dot{x}_{i}/2\!=\!\partial_{p_{i}}H$ and $\dot{p}_{i}/2\!=\!-\partial_{x_{i}}H$. Once $S(\mathbb{x})$ is known, by solving Hamilton’s equations, so is the MTE, which is proportional to $e^{NS(0,0)}$ KamenevPRE2008 ; Assaf2010 ; HindesPRL2016 .

For convenience, let us define new variables $u=(x_{1}-x_{2})/2$, $p_{u}=p_{1}-p_{2}$, $w=(x_{1}+x_{2})/2$ and $p_{w}=p_{1}+p_{2}$. This transformation is canonical since the determinant of the Jacobian $\partial(\textbf{Q},\textbf{P})/\partial(\textbf{x},\textbf{p})=1$, where $\textbf{Q}=(u,w)$, $\textbf{P}=(p_{u},p_{w})$, $\textbf{x}=(x_{1},x_{2})$, and $\textbf{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

[TABLE]

Since the transformation of variables is canonical, the action along the path to extinction is given by DykmanRev

[TABLE]

Transforming to the new variables in Hamiltonian (Degree Dispersion Increases the Rate of Rare Events in Population Networks), and assuming $u$ and $p_{u}$ scale as ${\cal O}(\epsilon)$, we find the trajectories $p_{w}(w)$ and $p_{u}(u)$ up to ${\cal O}(\epsilon^{2})$ SM . The trajectories are then substituted into Eq. (6), which yields

[TABLE]

where, $S_{0}=1/\Lambda+\ln\Lambda-1$ is the action for a degree-homogeneous network ($\epsilon=0$), and $f_{E}(\Lambda)>0$. 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 $\epsilon\ll 1$ our analytical results (Degree Dispersion Increases the Rate of Rare Events in Population Networks) agree well with numerical solutions of the Hamilton equations, obtained using the Iterative Action Minimization Method Lindley2 ; SM .

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 $\sum_{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}\!\equiv\!\langle k\rangle$, such that $g(k_{0}+i)=g(k_{0}-i)$ for $i=1,2,3,\dots$. Let us also assume our distribution has a bounded support such that $k_{\text{min}}=k_{0}-\Delta$ and $k_{\text{max}}=k_{0}+\Delta$, where $g(k<k_{\text{min}})=g(k>k_{\text{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_{\{n_{k}\}}$ – the joint probability to find $(n_{k_{\text{min}}},\dots,n_{k_{\text{max}}})$ infected nodes of degree $k$, and using the above WKB formalism, $P(\textbf{x})\sim e^{-NS(\textbf{x})}$, where $\textbf{x}=(x_{k_{\text{min}}},\dots,x_{k_{\text{max}}})$, we arrive at a Hamiltonian equivalent to HindesPRL2016 . Denoting $g(k)p_{k}\!=\!\partial S\!/\partial{x_{x}}$, the action can be shown to satisfy SM

[TABLE]

where we have used the symmetry of $g(k)$ about its mean $k_{0}$. Now, since each pair of nodes $k_{0}\pm j$ for $j\in[1,\Delta]$ can be viewed as a bimodal network, using Eqs. (6) and (Degree Dispersion Increases the Rate of Rare Events in Population Networks), the action for such a bimodal network with degrees $k_{0}-j$ and $k_{0}+j$, satisfies: $(1/2)\int p_{k_{0}-j}dx_{k_{0}-j}+p_{k_{0}+j}dx_{k_{0}+j}=S_{0}-f_{E}(\Lambda)\epsilon_{j}^{2}$, where $\epsilon_{j}\!=\!j/k_{0}$. Moreover, the node of rank $k_{0}$ can be viewed as a bimodal network with $\epsilon_{j}=0$, such that $\int p_{k_{0}}dx_{k_{0}}=S_{0}$. Therefore, using the fact that $\sum_{k}g(k)=1$ and that the variance of $g(k)$ satisfies $\sigma^{2}=\sum_{k}(k-k_{0})^{2}g(k)$, the action [Eq. (Degree Dispersion Increases the Rate of Rare Events in Population Networks)] and MTE become:

[TABLE]

Equation (9) is the first of the main results in this work. Namely for any network, if the CV is small, $\sigma/\langle k\rangle\!\ll\!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 $\sigma/\langle k\rangle\!\gg\!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 such as Salathe ; Vespignani1 . Furthermore, while the pre-factor for the relative increase of the logarithm of the MTE, $f_{E}(\Lambda)$, 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 $\sigma^{2}\!/\!\langle k\rangle^{2}$ and $\Lambda$.

Our analysis above required that the network degree distribution be symmetric and bounded. However, even for non-bounded 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 $\gamma_{1}$ 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 SM . 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 SM we show that our results even hold for power-law networks when the CV is not too largeSM .

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 Vespignani1 ; RednerBook . 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{\!\{\lambda\Delta E_{i}\}}]$, where $\Delta E_{i}$ is the change in the local pair-wise ferromagnetic energy for node $i$ to flip spin, and $\lambda$ is an inverse temperature. 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 $\mathbb{x}$ can be derived in precisely the same way as the SIS model above HindesSR2017 . The Hamiltonian reads

[TABLE]

where $\bar{x}\!=\!\sum_{k}\!kg(k)x_{k}/\!\left<k\right>$ is the degree-weighted mean magnetization, and $g(k)p_{k}\!=\!\partial S\!/\partial x_{k}$ are the momenta.

In contrast to the SIS model, the spin model exhibits three fixed points: $\mathbf{x}\!=\!\mathbf{x}^{*}$ and $\mathbf{x}\!=\!-\mathbf{x}^{*}$ which are stable, and $\mathbf{x}\!=\!\mathbf{0}$ which is unstable. The stable fixed points emerge at a pitchfork bifurcation when $\lambda=\lambda_{c}\!\equiv\!\left<k\right>\!/\!\left<k^{2}\right>$. As before, we may denote $\lambda\!=\!\Lambda\lambda_{c}$, where $\Lambda\!=\!1$ is the bifurcation threshold. In the spin model, demographic noise causes switching between $\mathbf{x}^{*}$ and $-\mathbf{x}^{*}$ ChenChaos . 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 FN_DB ). As a consequence, the deterministic trajectory starting from the vicinity of the unstable point $\mathbf{0}$ and ending at the stable fixed point $\mathbf{x}^{*}$, coincides up to time reversal, with the fluctuational path from $\mathbf{x}^{*}$ to $\mathbf{0}$ DykmanRev . Once at the unstable point $\mathbf{0}$, the network can switch to $-\mathbf{x}^{*}$ following its deterministic dynamics.

In order to find the switching path, we again use Hamilton’s equations $g(k)\dot{x}_{k}\!=\!\partial H/\partial p_{k}$. The relevant trajectories $p_{k}(\mathbf{x})$ can be found by equating $-\dot{x}_{k}|_{\mathbf{p}=\mathbf{0}}=\dot{x}_{k}(\mathbf{p})$, where the former represents the deterministic trajectory. By doing so, the switching path satisfies SM

[TABLE]

and hence the action for switching, $S(\mathbf{0})\!=\!\sum_{k}\!g(k)\int_{x_{k}^{*}}^{0_{k}}p_{k}dx_{k}$, becomes

[TABLE]

Following the same general approach as for the SIS model above, we write $k\!=\!k_{0}(1+\epsilon)$ where $\epsilon\equiv(k-k_{0})/k_{0}$. For degree distributions with a small CV, $\sigma/k_{0}\ll 1$, we have $\lambda\!\approx\!\Lambda[1-\left<\epsilon^{2}\right>]/k_{0}$ and $\left<\epsilon^{2}\right>=\sigma^{2}/k_{0}^{2}$, as before. In order to evaluate Eq. (11) in the limit of $\left<|\epsilon|\right>\!\ll\!1$, we use the small-$\left<|\epsilon|\right>$ expansion of $x_{k}^{*}$ and $\bar{x}^{*}$, see SM , and keep terms up to order $\left<\epsilon^{2}\right>$. This procedure yields the action and mean switching time (MST)

[TABLE]

where $S_{0}=-(1/2)\left[\ln\left(1-x_{0}^{2}\right)+\Lambda x_{0}^{2}\right]>0$, $x_{0}$ is the positive solution of $x_{0}\!=\!\tanh\{\Lambda x_{0}\}$, and $f_{S}(\Lambda)>0$.

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. (Degree Dispersion Increases the Rate of Rare Events in Population Networks) 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(\mathbf{0})-S_{0}]/f(\Lambda)$ versus the CV, and obtain a collapse across all networks and all $\Lambda$, for both models: network simulations and numerical solutions of the Hamilton equations SM . As our analysis exemplifies, if the rate of rare events (on log scale) is normalized by the correct process-dependent factor, $f(\Lambda)$, 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 SM for power-law networks and continuous-noise analogs for both processes SM .

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 sourcesSM . 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 bifurcationSM .

We thank Lev Muchnik and Ira B. Schwartz for useful discussions, and Baruch Meerson for critically reading the manuscript. MA was supported through the Israel Science Foundation Grant No. 300/14 and the United States-Israel Binational Science Foundation grant No. 2016-655. JH was supported through the U.S Naval Research Laboratory Karle Fellowship.

Supplemental Material

I SIS Hamiltonian for arbitrary degree distributions

Following the main text, we first write a master equation for ${\bf n}$, where $n_{k}$ is the number of infected nodes with degree $k$, $N$ is the total number of nodes in the network, $N_{k}$ is the total number of nodes of degree $k$, and $g(k)=N_{k}/N$ is the node degree distribution. Given the annealed network approximation and current state ${\bf n}$, the rate at which $n_{k}$ increases by one is $\lambda k(N_{k}-n_{k})\bar{x}$, where $\bar{x}\equiv\sum_{k^{\prime}}k^{\prime}n_{k}^{\prime}/[N\!\left<k\right>]=\sum_{k^{\prime}}g_{k^{\prime}}k^{\prime}x_{k^{\prime}}/\!\left<k\right>$, is the fraction of infected neighbors along an edge, and $x_{k}=n_{k}/N_{k}$ is the fraction of infected nodes of degree k. Similarly, the rate at which $n_{k}$ is decreased by one is $n_{k}$. We can denote these transitions, compactly, with the notation ${\bf n}\rightarrow{\bf n}+\bm{1}_{k}$ and ${\bf n}\rightarrow{\bf n}-\bm{1}_{k}$, respectively.

Consequently, the master equation reads

[TABLE]

Now we assume the system has entered a long-lived metastable state, such that $\partial P({\bf n},t)/\partial t\simeq 0$, use the WKB ansatz for the quasi-stationary distribution $P({\bf n})\sim e^{-NS({\bf x})}$, and keep only leading-order terms in $N\gg 1$. This gives rise to a stationary Hamilton-Jacobi equation (where the action has no explicit time dependence), $H({\bf x},\partial_{{\bf x}}S)=0$, with

[TABLE]

The momenta, $\lambda_{k}\equiv\partial_{{x}_{k}}S$, can be usefully redefined as $p_{k}=\lambda_{k}/g_{k}$. With this transformation $\dot{x}_{k}=\partial_{\lambda_{k}}H({\bf x},{\bf\lambda})=\partial_{p_{k}}H({\bf x},{\bf p})/g_{k}$. Similarly, since $\dot{\lambda}_{k}=-\partial_{x_{k}}H({\bf x},{\bf\lambda})$, we get $\dot{p}_{k}=-\partial_{x_{k}}H({\bf x},{\bf p})/g_{k}$. As a result, the action satisfies

[TABLE]

We note that in Ref. [21] “$y_{k}$” is what we call $x_{k}$ in this work.

II Finding the optimal path in the SIS model

In this section we consider a bimodal network with only two degrees $k_{1}=k_{0}(1-\epsilon)$ and $k_{2}=k_{0}(1+\epsilon)$, where $\langle k\rangle=k_{0}$ is the mean degree of the network, $\sigma=\epsilon k_{0}$ is its standard deviation, while $\epsilon\ll 1$. Following the main text, here we find the optimal path to extinction, and the action along it, for such a bimodal network.

To conveniently deal with the Hamiltonian [Eq. (4) in the main text] in the limit $\epsilon\ll 1$, let us define new variables $u=(x_{1}-x_{2})/2$, $p_{u}=p_{1}-p_{2}$, $w=(x_{1}+x_{2})/2$ and $p_{w}=p_{1}+p_{2}$. This transformation is canonical since the determinant of the Jacobian $\partial(\textbf{Q},\textbf{P})/\partial(\textbf{x},\textbf{p})=1$, where $\textbf{Q}=(u,w)$, $\textbf{P}=(p_{u},p_{w})$, $\textbf{x}=(x_{1},x_{2})$, and $\textbf{p}=(p_{1},p_{2})$. Using the new variables, the path to extinction is a heteroclinic trajectory (or instanton) connecting between the fixed points $[w^{*},u^{*},0,0]$ and $[0,0,p_{w}^{*},p_{u}^{*}]$, where

[TABLE]

and $x_{0}=(\Lambda-1)/\Lambda$. Since the transformation of variables is canonical, the action along the path to extinction is given by [1]

[TABLE]

In the following we find the trajectories $p_{w}(w)$ and $p_{u}(u)$, and compute the integral (17). We begin by finding $p_{w}(w)$. Plugging $x_{1}=w+u$, $x_{2}=w-u$, $p_{1}=(p_{w}+p_{u})/2$ and $p_{2}=(p_{w}-p_{u})/2$ into the Hamiltonian [Eq. (4) in the main text], and assuming $u$ and $p_{u}$ scale as ${\cal O}(\epsilon)$, we find in the leading order

[TABLE]

As a result, we find in the leading order $p_{w}^{(0)}=-2\ln[\Lambda(1-w)]$. To find the subleading ${\cal O}(\epsilon^{2})$ correction, we demand that (i) $p_{w}$ vanish at $w=w^{*}$, and (ii) $p_{w}=p_{w}^{*}$ at $w=0$. If we simply interpolate between the two fixed points of $p_{w}(w)$ by using a linear function of $w$, we get

[TABLE]

One can check a-posteriori that $p_{w}(0)=p_{w}^{*}$ and $p_{w}(w^{*})=0$ up to ${\cal O}(\epsilon^{4})$ corrections. In Fig. 5 we numerically verify that Eq. (19) holds up to ${\cal O}(\epsilon^{2})$. Note, that the numerical solutions of the Hamilton equations, which yield the optimal paths to extinction/switching and the corresponding actions along these paths, were found by using the Iterative Action Minimization Method, see Ref. [26] for further details. Matlab code is available upon request.

Regarding $p_{u}(u)$, we notice that both $p_{u}^{*}$ and $u^{*}$ scale as ${\cal O}(\epsilon)$, and thus we expect both $u$ and $p_{u}$ to scale as ${\cal O}(\epsilon)$ in the entire path. Since the integral over $p_{u}du$ already scales as ${\cal O}(\epsilon^{2})$, it is sufficient to approximate $p_{u}(u)$ as a straight line connecting $[u^{*},0]$ and $[0,p_{u}^{*}]$:

[TABLE]

which vanishes at $u=u^{*}$ and equals $p_{u}^{*}$ at $u=0$. Again, this choice of path agrees well with numerics, see Fig. 5.

Finally, performing the integrations in Eq. (17) using Eqs. (19) and (20) and keeping terms up to ${\cal O}(\epsilon^{2})$, gives

[TABLE]

where $S_{0}=1/\Lambda+\ln\Lambda-1$ is the action for a degree-homogeneous network ($\epsilon=0$).

III Extension of the SIS result to non-symmetric distributions

Here we generalize Eq. (9) in the main text to non-symmetric degree distributions. For any degree distribution, the action along the optimal path is given by

[TABLE]

Let us assume a general distribution $g(k)$ centered about $k_{0}$, with $\sigma\ll k_{0}$. Thus, it is sufficient to take the sum up to $2k_{0}$, since the width is much smaller than the mean and $g(2k_{0})$ is already negligible. Denoting by $I_{k}=\int p_{k}dx_{k}$, we have

[TABLE]

where $\delta_{j}=[g(k_{0}+j)-g(k_{0}-j)]/g(k_{0}-j)$ denotes the deviation from symmetry of the degree distribution. Taylor-expanding $\delta_{j}$ around $k_{0}$ up to third order, we find $\delta_{j}\simeq[g^{\prime\prime\prime}(k_{0})/(3g(k_{0}))]j^{3}$, where prime denotes differentiation with respect to the degree $k$. Evaluating this term at $j=\sigma$, where the distribution has already decayed by a factor of $e$, we find $\delta_{j}(j=\sigma)\simeq[g^{\prime\prime\prime}(k_{0})/(3g(k_{0}))]\sigma^{3}$. We have evaluated this term for various examples of degree distributions including the Poisson and Gamma distributions, and found in all examples that $\delta_{j}(j=\sigma)$ is proportional to the distribution’s skewness $\gamma_{1}$. Therefore, for distributions with a small skewness, $\delta_{j}\sim\gamma_{1}\ll 1$ for $j={\cal O}(\sigma)$. For smaller $j$, obviously $\delta_{j}$ is smaller (and again negligible compared to $1$), as we are in the symmetric region of the distribution, while for $j\gg\sigma$, the distribution has already decayed and the terms in the sum are negligible. As a result, we can safely neglect $\delta_{j}$ in Eq. (III) for all $j$’s, and we recover Eqs. (8) and (9) in the main text, which were derived for symmetrical distributions.

IV Finding the optimal path and action in the spin model

Here we consider the spin model and find the switching path (or instanton) along which the action can be calculated. To do so, we use Hamilton’s equations $g(k)\dot{x}_{k}\!=\!\partial H/\partial p_{k}$, where the Hamiltonian is given by Eq. (10) in the main text. The relevant trajectories $p_{k}(\mathbf{x})$ can be found by equating $-\dot{x}_{k}|_{\mathbf{p}=\mathbf{0}}=\dot{x}_{k}(\mathbf{p})$, where the former represents the deterministic trajectory. By doing so, we obtain the following equations for $p_{k}(\mathbf{x})$:

[TABLE]

After some algebra, we find a solution

[TABLE]

which leads to the action [Eq. (11) in the main text].

In order to approximate the action in the limit of $\left<|\epsilon|\right>=\sigma/k_{0}\!\ll\!1$ we need to first evaluate $x_{k}^{*}$ and $\bar{x}^{*}$ in that limit. Using the Hamiltonian [Eq. (10) in the main text], the deterministic rate equations (when $\mathbf{p}\!=\!\mathbf{0}$) have fixed points $x_{k}^{*}$ which satisfy the following transcendental equations: $x_{k}^{*}\!=\!\tanh\{\lambda k\bar{x}^{*}\}$ [36,37]. If we assume that $\bar{x}^{*}$ takes the form $\bar{x}^{*}\!=\!x_{0}+\left<\epsilon^{2}\right>x_{1}$, where $x_{0}$ is the positive solution of $x_{0}\!=\!\tanh\{\Lambda x_{0}\}$, then

[TABLE]

Substituting Eq. (26) into the definition of $\bar{x}^{*}$ we find

[TABLE]

where $x_{1}$ in Eq. (26) satisfies $x_{1}=-\Lambda^{2}x_{0}^{3}(1-x_{0}^{2})/[1-(1-x_{0}^{2})\Lambda]$. Plugging Eqs. (26) and (27) into the action [Eq. (11) in the main text] yields the final result for the mean switching time

[TABLE]

where $S_{0}=-(1/2)\left[\ln\left(1-x_{0}^{2}\right)+\Lambda x_{0}^{2}\right]>0$. This result coincides with Eq. (12) in the main text.

V Breaking detailed balance in the spin model

Here we generalize our results for the spin model in the absence of detailed balance. A simple way to break detailed balance is to add a spontaneous transition with rate $f$. Namely, we assume that each node flips spin at a stochastic rate, $f+[1+\exp{\!\{\lambda\Delta E_{i}\}}]^{-1}$. In the presence of this spontaneous flipping process, the Hamiltonian [Eq. (10) in the main text] becomes:

[TABLE]

which can be derived in exactly the same way as above for the SIS model (see Ref. [37]). It is straightforward to show that the pitchfork bifurcation now occurs when $\lambda\left<k^{2}\right>\!/\!\left<k\right>-1-2f\!=\!0$.

The action for switching can be computed from

[TABLE]

where $g_{k}\dot{x}_{k}\!=\!\partial H/\partial p_{k}$ and $g_{k}\dot{p}_{k}\!=-\partial H/\partial x_{k}$. We solve this system numerically for several networks and values of $f$; the results are shown in Fig. 6. In order to keep the distance to bifurcation constant across all networks used, we define $\lambda\!=\![1+\delta]\left<k\right>\!/\!\left<k^{2}\right>$. Therefore, all three series in Fig. 6 have the same distance to bifurcation, $\lambda\left<k^{2}\right>\!/\!\left<k\right>-1-2f$.

Our numerical results indicate that, even in the absence of detailed balance, the correction to the action across all networks collapses to the same expression

[TABLE]

and hence, our main result is preserved. Note however, that $f_{S}(\delta,f)$ is no longer a function, only, of the distance to bifurcation – otherwise all three series would collapse to the same correction. This more general function could be calculated, i.e., with the general procedure used for extinction in the SIS model, without assuming detailed balance; see main text and Sec. SM-II.

VI Parameters for Fig. 4 in main text

Here we describe in detail the results shown in Fig. 4 in the main text. The network simulations for this figure were taken from Fig. 2 (left) and Fig. 3 (left) in the main text. Furthermore, Fig. 4 in the main text includes numerical solutions of Hamilton’s equations. In red we show the numerical results for extinction; circles are bimodal distributions with $\left<k\right>\!=\!40$ and $\Lambda\!=\!3.5$, squares are uniform distributions with $\left<k\right>\!=\!30$ and $\Lambda\!=\!2.0$, crosses are generalized Gaussian distributions with exponent $1.0$, $\left<k\right>\!=\!35$, and $\Lambda\!=\!3.0$. In blue, we show numerical results for switching; circles are bimodal distributions with $\left<k\right>\!=\!50$ and $\Lambda\!=\!1.4$, squares are uniform distributions with $\left<k\right>\!=\!50$ and $\Lambda\!=\!1.4$, diamonds are Gamma distributions with $\left<k\right>\!=\!50$ and $\Lambda\!=\!1.4$, triangles are Gaussian distributions with $\left<k\right>\!=\!108.5$ and $\Lambda\!=\!1.7$, and crosses are generalized Gaussian distributions with exponent $0.75$, $\left<k\right>\!=\!35$, and $\Lambda\!=\!1.7$. Note that given these parameters, each degree distribution has a single parameter which can be varied to change the coefficient of variation. Finally, the dashed-line in Fig. 4 in the main text is the theoretical prediction $-[S({\bf 0})-S_{0}]/f(\Lambda)=\sigma^{2}/\langle k\rangle^{2}$.

 

Here and in Figs. 2 and 3 in the main text the simulations on networks were performed using Monte-Carlo simulations implemented according to the Gillepsie’s algorithm in continuous time. Namely, for each node in a network, there is an exponentially distributed time to make a transition to another state. For example, a transition of a susceptible node to infected occurs at a rate $\lambda$ times the number of infected neighbors. Noise comes from the fact that the time is not deterministic, but is a stochastic variable. C++ code is available upon request.

VII Continuous models with continuous noise

In the main text, we deal with discrete states on the nodes. However, our results are qualitatively the same for continuous states with continuous noise, and an analogous perturbation-theory in $\sigma\!/\!\left<k\right>$ can be developed. In particular the network action for extinction/switching takes the form, S\big{(}\Lambda,\frac{\sigma}{\left<k\right>}\big{)}\approx S\big{(}\Lambda,0\big{)}-f(\Lambda)\frac{\sigma^{2}}{\left<k\right>^{2}}.

Let us consider the following Langevin system

[TABLE]

where $\mathcal{M}_{i}(\mathbf{x};\lambda)$ is the mean-field dynamics for node $i$, and $\xi_{i}(t)$ is independent and identically distributed Gaussian white noise (GWN), $\left<\xi_{i}(t)\xi_{j}(t^{\prime})\right>=2D\delta_{ij}\delta(t-t^{\prime}).$ The mean-field dynamics correspond to $\mathbf{p}\rightarrow\mathbf{0}$ in Hamilton’s equations, or

[TABLE]

for the SIS and spin models, respectively.

Similar to the main text, the quasi-stationary probability distributions have a WKB form when $D\ll 1$, $P(\mathbf{x})\!\sim\!\exp\{-NS(\mathbf{x})/[2D]\}$. Hamilton’s equations are straightforward to derive (see for instance E. Forgoston and R. O. Moore, SIAM Rev. 60(4), 969 (2018)), and represent an application of classical large-deviation theory for dynamical systems perturbed by GWN. The procedure for deriving Hamilton’s equations is essentially the same as in the main-text and Sec. SM-I, except the master equation, e.g. Eq.(I), is replaced by a Fokker-Planck equation for Eq.(32).

Given that we expect nodes with the same degree to have synchronized dynamics during a large fluctuation (i.e., trading the node subscript $i$ for the degree subscript $k$), we find for the SIS model

[TABLE]

and for the spin model

[TABLE]

where, as above, $\bar{x}\!\equiv\!\sum_{k^{\prime}}\!k^{\prime}x_{k^{\prime}}/[N\!\left<k\right>]$ and $S\!=\!\sum_{k}g_{k}\int p_{k}dx_{k}$.

Figure 7 shows the change in the action from the homogeneous network limit, for both processes, as a function of $\sigma^{2}\!/\!\left<k\right>^{2}$ for three different examples of degree distributions. The results are consistent with those presented in the main text.

VIII Power-law degree distributions

In the main-text we primarily discuss networks whose degree-distributions are centralized around a mean, with an approximately symmetric pattern of dispersion. Nevertheless, our quantitative results turn out to also hold for power-law networks with relatively small coefficients of variation (e.g. degree exponents greater than four). Moreover, our qualitative result: degree dispersion increases the rate of rare events when comparing networks with constant distances to threshold, holds for power-law networks with even smaller degree exponents. Figure 8 shows the action for extinction for power-law networks with degree distributions $g(k,s)\!=\!k^{-s}\!/\!\sum_{k^{\prime}=20}^{500}k^{\prime-s}$. The degree exponent ranges from $s\!=\!10,9.5,...,2.5$, where $\Lambda\!=\!1.5$ (red) and $\Lambda\!=\!2.0$ (blue). The dashed line shows the predicted scaling, $S(\mathbf{0},s)\!\approx\!S(\mathbf{0},s\!\rightarrow\!\infty)-f_{E}(\Lambda)\frac{\sigma^{2}}{\left<k\right>^{2}}$ (Eq.(7) and Eq.(9) in the main text), which agrees well with numerics for $s\!>\!5$. For reference, a power-law network with $s\!=\!5$ has a variance of $\sigma^{2}\!/\!\left<k\right>^{2}\!\approx\!0.12$.

IX Generality of our results

In this section we briefly discuss the generality of our results. We have shown that the barrier for extinction/switching, given in the form of a cumulative action obtained by integrating over a trajectory between the deterministically stable and unstable fixed points, decreases as the heterogeneity of the network, $\epsilon\!=\!\sigma\!/\!\left<k\right>$, is increased. Specifically, we have demonstrated the following functional dependence S\big{(}\Lambda,\epsilon\big{)}\approx S\big{(}\Lambda,0\big{)}-f(\Lambda)\epsilon^{2}, for both the SIS model of epidemics and a model of spontaneous magnetization flipping, where $f(\Lambda)$ depends on the local microscopic dynamics, but is independent on the network topology. That is, as long as the heterogeneity parameter $\epsilon$ is fixed, we have shown that the network topology affects the mean escape time in a universal manner, regardless of the degree distribution of the network. Moreover, the dependence on heterogeneity holds for distinct types of rare events: extinction and switching.

It is our conjecture, that any model that satisfies the following generic conditions will demonstrate similar quantitative features:

• •

At the microscopic level, the model should include one-body and two-body interactions, where the latter are due to interactions between each node and its neighbors. These microscopic dynamics determine the specific nature of the function $f(\Lambda)$.

• •

The microscopic dynamics should give rise at the deterministic level to a nontrivial stable state and an adjacent unstable state which is either an absorbing state, or it is accompanied by an additional (target) stable fixed point. Such states should be fixed-points of a mean-field dynamics. The mean-field description should entail sets of differential equations in time for the density, or set of densities, describing the average state of nodes with degree $k$ (or eigenvector centrality).

• •

The model should include a tuning parameter $\Lambda$, which when approaches some $\Lambda_{c}$, the stable fixed point(s) at the deterministic level merge with the unstable fixed point, and the system becomes monostable. That is, the system can undergo, e.g., a transcritical, a pitchfork or a saddle-node bifurcation, depending on the scenario at hand. When noise is accounted for, the former case typically corresponds to an escape from a metastable state to an absorbing state (e.g., extinction), while the latter cases typically correspond to switching between two metastable states separated by a saddle point.

Note: the state space and noise can be continuous or discrete.

Finally, while we have considered two prototypical examples of extinction and switching, we expect our results to hold for wide variety of additional models which satisfy the conditions specified above. Examples include population dynamics models (or equivalent) with an Allee effect, other models of epidemics such as the SIRS models, and voter models on networks with hysteresis. On the other hand, models of evolutionary game theory on networks and generalized contagion models, which include more complicated bifurcation scenarios, are expected to (possibly) display a different dependence on $\epsilon$ in the action, as the network heterogeneity is increased.