Epidemic Extinction Paths in Complex Networks

TL;DR

This paper analyzes the stochastic extinction of epidemics in complex networks, predicting the most probable paths to disease eradication and quantifying fluctuation patterns using analytical and simulation methods.

Contribution

It introduces an analytical framework for predicting epidemic extinction paths and times in complex networks, validated by simulations across various network types.

Findings

Predicted distribution of large fluctuations leading to extinction.

Identified most probable extinction paths in different network structures.

Quantified scaling patterns of extinction phenomena near thresholds.

Abstract

We study the extinction of long-lived epidemics on finite complex networks induced by intrinsic noise. Applying analytical techniques to the stochastic Susceptible-Infected-Susceptible model, we predict the distribution of large fluctuations, the most probable, or optimal path through a network that leads to a disease-free state from an endemic state, and the average extinction time in general configurations. Our predictions agree with Monte-Carlo simulations on several networks, including synthetic weighted and degree-distributed networks with degree correlations, and an empirical high school contact network. In addition, our approach quantifies characteristic scaling patterns for the optimal path and distribution of large fluctuations, both near and away from the epidemic threshold, in networks with heterogeneous eigenvector centrality and degree distributions.

Figures

Equations

$$\frac{\partial P}{\partial t}(\mathbb{I},t)=\sum _i{R}_i^{+}(\mathbb{I}-1_i)P(\mathbb{I}-1_i,t)-R_i^{+}(\mathbb{I})P(\mathbb{I},t)$$

$$+R_i^{-}(\mathbb{I}+1_i)P(\mathbb{I}+1_i,t)-R_i^{-}(\mathbb{I})P(\mathbb{I},t),$$

$$\rho (\mathbf{\nu },t)\cong \rho (x,t)=b{e}^{-S(x,t)}.$$

$$\partial S/\partial t+H(x,\partial S/\partial x)=0,$$

$$\displaystyle\!\!H(\mathbb{x},\mathbb{p})\!=\!\sum_{i}\!\!\Bigg{[}\!\beta\big{(}1\!-\!x_{i}\big{)}\!\big{(}e^{p_{i}}\!-\!1\big{)}\!\!\sum_{j}\!A_{ij}x_{j}+\alpha x_{i}\big{(}e^{-p_{i}}\!-\!1\big{)}\!\Bigg{]}\!.$$

$${\dot{x}}_i$$

$${\dot{p}}_i$$

$$S(x,t)={\int }_{x(t=t_0)}^{x}p\cdot dx-{\int }_{t_0}^{t}H(x,p)d{t}^{\prime },$$

$$⟨T⟩=B(\tilde{\beta },A)e^{S(x=0)}/\alpha ,$$

$${\lambda }^{(1)}{x}_{i,1}$$

$${\lambda }^{(1)}{x}_{i,2}$$

$$x_i^{\ast }$$

$$p_i^{\ast }$$

$$S(x(a))=\frac{{\delta }^2\left(a-a^2/2\right)}{{\sum }_j{\zeta }_j^{(1)}{{\eta }_j^{(1)}}^2{\sum }_l{\eta }_l^{(1)}{{\zeta }_l^{(1)}}^2}+\mathcal{O}({\delta }^3).$$

$${\dot{x}}_b=$$

$${\dot{p}}_b=$$

$$-e^{-p_b}+1.$$

$$x_i^{\ast }$$

$$p_i^{\ast }$$

$$\displaystyle 1=\sum_{j}\tilde{\beta}\lambda\eta_{j}^{2}\big{(}1-x^{*}_{j}\big{)}\big{/}\big{[}1+\tilde{\beta}\lambda\eta_{j}X-\sigma^{o}\big{]},$$

$$\displaystyle\mu_{i}^{o}=\tilde{\beta}\lambda\eta_{i}\!\sum_{j}\eta_{j}\mu_{j}^{o}(1-x_{j}^{*})/\big{[}1+\tilde{\beta}\lambda\eta_{i}X-\sigma^{o}\big{]}.$$

$$\displaystyle 1=\sum_{j}\tilde{\beta}\lambda\eta_{j}^{2}e^{p_{j}^{*}}\big{/}\big{[}\sigma^{in}+e^{-p_{j}^{*}}\big{]},$$

$$\displaystyle\epsilon_{i}^{in}=\tilde{\beta}\lambda\eta_{i}e^{p_{i}^{*}}\!\sum_{j}\eta_{j}\epsilon_{j}^{in}/\big{[}\sigma^{in}+e^{-p_{i}^{*}}\big{]}.$$

$${\dot{x}}_k=$$

$${\dot{p}}_k=$$

$$-e^{-p_k}+1.$$

$$S(x)=N\sum _{k}g(k){\int }_{x_k^{\ast }}^{x_k}{p}_{k}d{x}_k^{\prime }.$$

$$S=\frac{1}{2}N{\delta }^2\frac{{⟨v^2⟩}^3}{{⟨v^3⟩}^2}+\mathcal{O}({\delta }^3),$$

$$S=\frac{{\delta }^2}{2}\frac{(1-\alpha {)}^2(2+\alpha )}{{\alpha }^3}{N}^{1-2[1-\alpha ]/[2+\alpha ]},$$

$$\displaystyle\sum_{j}\frac{\eta_{j}\epsilon_{j}^{o}}{N\!\left<\eta\right>}\approx\;\epsilon_{i}^{o}\Big{[}\sigma^{o}\!+\!1\!+\tilde{\beta}\lambda\eta_{i}N\!\left<\eta\right>\!\Big{(}1\!-\!\big{(}1\big{/}\tilde{\beta}\lambda N\!\left<\eta\right>^{\!2}\!\big{)}\!\Big{)}\!\Big{]}$$

$$\displaystyle-\mu_{i}^{o}\Big{[}2\!-\!\big{(}1\big{/}\tilde{\beta}\lambda N\!\left<\eta\right>^{\!2}\!\big{)}\!-\!\big{(}1\big{/}\tilde{\beta}\lambda\eta_{i}N\!\left<\eta\right>\!\big{)}\Big{]}.$$

$$\displaystyle\left[-1+\big{(}(\sigma^{in}-1)\!\big{/}\tilde{\beta}\lambda\eta_{i}N\!\left<\eta\right>\!\big{)}\!\right]\!\mu_{i}^{in}\approx\;\!\!-\!\sum_{j}\!\frac{\eta_{j}\mu_{j}^{in}}{N\!\left<\eta\right>}\frac{1}{\tilde{\beta}\lambda\eta_{j}N\!\left<\eta\right>}$$

$$\displaystyle+\!\sum_{j}\frac{\eta_{j}\epsilon_{j}^{in}}{N\!\left<\eta\right>}\!\left[-2+(1\big{/}\tilde{\beta}\lambda\eta_{i}N\!\left<\eta\right>)+(1\big{/}\tilde{\beta}\lambda\eta_{j}N\!\left<\eta\right>)\right].$$

$${\sigma }_i^{in}$$

$${\sigma }_i^o$$

Full text

I INTRODUCTION

Understanding the dynamics of infectious processes in complex networks is an important problem, both in terms of generalizing concepts in statistical mechanics and applying them to public health Pastor ; Anderson ; Dorogovtsev ; Hindes1 ; Vespignani1 . A primary question in infectious disease modeling is how to control an outbreak, with the ultimate goal of reducing the number of individuals able to spread infection to zero. This process, by which an epidemic is extinguished, is called extinction or disease fade-out Doering ; Schwartz ; Meerson2 ; Meerson3 ; Assaf2 . To understand and possibly achieve extinction, mathematical models can be useful, where extinction can be naturally captured in terms of a dynamical transition from an endemic state (e.g., fluctuating equilibrium or cycle) to a disease-free state Assaf1 .

Although it is known that random fluctuations are the cause of extinction in finite populations, the process of extinction does not happen in the deterministic systems analyzed in the vast majority of works on endemic dynamics in networks – where contacts between infectious and susceptible individuals are typically assumed to be well above an epidemic threshold or bifurcation point Pastor ; Meerson3 ; Nasell . Consequently network-control prescriptions often reduce to bringing systems below a bifurcation point VaccRev ; Schwartz ; Drakopoulos . One may ask, is targeting sub-threshold regimes as a control method necessary or even optimal? In actuality the spread of disease is a highly stochastic process both in terms of the natural randomness inherent in contact processes and fluctuations due to time-varying and uncertain environments Meerson3 ; Meerson . These stochastic effects make extinction inevitable, even above threshold, in finite networks and should be reflected in epidemic controls Schwartz2 ; Billings . In fact, recent work has shown that optimal control of networks with noisy dynamics leverages randomness and a network’s natural, noise-induced pathway between distinct states Hindes1 ; Motter . Continuing in this line of thinking, we seek a prescription for computing epidemic extinction pathways through complex networks.

Such issues have received much attention in well-mixed and spatially homogeneous models Dykman ; Meerson ; Kamenev ; Billings ; Schwartz2 ; Meerson4 . It has been demonstrated in many works that noise and a system’s dynamics can couple in such a way as to induce a large fluctuation – effectively driving a system from one state to another Schwartz3 . If the fluctuation is a rare event in the weak noise limit, then the process is captured by a path that is a maximum in probability, or optimal path (OP), where all others are exponentially less likely to occur. The formalism borrows from analytical mechanics, describing the OP as a least-action trajectory in some effectively classical system, and allows one to predict the dynamical extinction pathway and the average time needed to realize it Dykman2 ; Friedlin . Some recent works have made progress in understanding extinction in networks, (e.g., deriving bounds for average extinction times), but do not make use of the path-based formalism outlined here Durrett ; Mountford ; Van ; Munoz ; Buono ; Assaf3 .

The following layout of the paper describes epidemic extinction through complex networks with intrinsic (demographic) noise in terms of large fluctuations and rare-event theory. Sec.II constructs the formalism: combining a mean-field approximation for endemic dynamics on networks with a Wentzel–Kramers–Brillouin $(WKB)$ technique that allows for an analytical description of the distribution of large fluctuations and the OP. The limiting form of the OP is discussed near the epidemic threshold in Sec.II.1, and away from threshold in Sec.III.1.1. Sec.III addresses how to compute the OP, extinction time, and their dependencies in certain cases, including in networks with large spectral gaps (Sec.III.1) and degree distributions (Sec.III.2). Throughout, predictions are compared to real and synthetic network simulations.

II LARGE FLUCTUATIONS, MEAN-FIELD, AND WKB APPROXIMATION

In order to predict epidemic extinction in general contact networks it is necessary to consider an arbitrary weighted adjacency matrix, $A$, where each element, $A_{ij}$, represents the strength of a link, or contact, from node $i$ to node $j$, in a graph with $N$ nodes. Given this representation, a network’s epidemic dynamics, assuming a simple Susceptible-Infectious-Susceptible Markov process (SIS) is captured by the states and transitions of its nodes; i.e., node $i$ is either “infected”, denoted $\nu_{i}\!=\!1$, or “susceptible”, $\nu_{i}\!=\!0$. Furthermore, node $i$ changes its state $\nu_{i}\!:\!0\rightarrow\!1$ with probability per unit time $\beta(1-\nu_{i})\sum_{j}A_{ij}\nu_{j}$, and $\nu_{i}\!:\!1\rightarrow\!0$ with probability per unit time $\alpha\nu_{i}$, where $\beta$ and $\alpha$ are known as the infection and recovery rates, respectively Pastor ; Hindes1 ; Vespignani ; Vespignani1 . Since the elements of $A$ are proportional to probabilities, $A$ is assumed to be nonnegative. It is important to note that there is inherent noise in the SIS model defined, which arises from the underlying stochastic reactions Gillespie .

In order to analyze the stochastic dynamics, it is useful to consider an ensemble consisting of $C$ identical networks with the same $A$, but independent realizations of the stochastic dynamics Ott2 . Each node can be specified by a graph position $i$, and ensemble number $c$, with state $\nu_{i,c}\!\in\!\{0,1\}$. In this way, the number of infected nodes in the ensemble with graph position $i$ is $I_{i}\!=\!\sum_{c}\nu_{i,c}$, with corresponding transitions and rates: $I_{i}\!\rightarrow\!I_{i}+1$ with rate $R_{i}^{+}(\mathbb{I})\!\equiv\!\beta\sum_{c}(1-\nu_{i,c})\sum_{j}A_{ij}\nu_{j,c}\!=\!\beta\sum_{j}A_{ij}\sum_{c}(1-\nu_{i,c})\nu_{j,c}$, and $I_{i}\!\rightarrow\!I_{i}-1$ with rate $R_{i}^{-}(\mathbb{I})\!\equiv\!\alpha I_{i}$. To simplify our analysis, it is useful to make a mean-field approximation and replace $\nu_{i,c}$ by the ensemble average, $I_{i}/C$: $R_{i}^{+}(\mathbb{I})\!\approx\!\beta\sum_{j}A_{ij}I_{j}(1-I_{i}/C)$, so that the transition rates depend explicitly on $\mathbb{I}=\left<I_{1},I_{2},...,I_{N}\right>$ alone. This approximation neglects correlations between neighboring graph positions Goltsev ; Mata ; Pastor ; FN5 . Ultimately, we are interested in the limit of large $C$, so that $x_{i}\!\equiv\!I_{i}/C$ gives a continuous fraction of infected nodes, or density, in graph position $i$. In this way, the large ensemble allows us to consider continuous densities even in discrete networks with unique graph positions.

Given the stochastic reactions and rates $R_{i}^{+}$ and $R_{i}^{-}$, the ensemble dynamics is described by a probability distribution, $P(\mathbb{I},t)$, satisfying a master equation:

[TABLE]

where $\mathbb{1}_{i}=\left<0\;{}_{1},0\;{}_{2},...,1\;{}_{i},0\;{}_{i+1},...\right>$. Because extinctions in large networks ($N\!\gg\!1$) with long-lived epidemics are rare events with small probabilities, we are interested in the tails of $P(\mathbb{I},t)$, where $\mathbb{I}$ corresponds to a large deviation from the average behavior, and is accompanied by an exponential reduction in probability.

This intuition suggests looking for solutions of Eq.(1) with an exponential, or $WKB$ form, $P(\mathbb{I},t)=ae^{-CS(\mathbb{x},t)}$ Assaf1 ; Dykman2 . The WKB solution for the ensemble distribution can be viewed as a product of independent and identical distributions for each realization in the ensemble. Hence, we can approximate the probability distribution of states for a single realization, $\rho(\bm{\nu},t)$, by

[TABLE]

Predictions from Eq.(2) (when combined with Eqs.(5-7) below) are in good agreement with simulations on an empirical high school network Salathe , shown in red in Fig.1.

We can find the leading contribution to $P(\mathbb{I},t)$ by substituting the WKB ansatz into Eq.(1), expanding in powers of the small parameter $1/C$ (e.g., $S(\mathbb{x}\pm\mathbb{1}_{i}/C)\!\approx\!S(\mathbb{x})\pm(1/C)\partial S/\partial x_{i}$), and neglecting terms of $\mathcal{O}(1/C)$ or smaller, where $C\!\gg\!1$. This approximation converts the master equation into a Hamilton-Jacobi equation (HJE):

[TABLE]

where $S$ and $H$ are called the Action and Hamiltonian, respectively. The latter is a function of the infection density at graph position $i$, $x_{i}$, and its conjugate momentum, $p_{i}=\partial S/\partial x_{i}$:

[TABLE]

Just as in analytical mechanics, a convenient approach for solving the HJE is to solve Hamilton’s equations of motion, $\dot{x}_{i}=\partial H/\partial p_{i}$ and $\dot{p}_{i}=-\partial H/\partial x_{i}$:

[TABLE]

expressed in terms of the ratio $\tilde{\beta}\!\equiv\!\beta/\alpha$, and the time, $\tau\!\equiv\!\alpha t$. Crucially, solutions of the HJE extremize $S$, when expressed as the integral

[TABLE]

where $\mathbb{x}(t)$ and $\mathbb{p}(t)$ are determined from Eqs.(5-6)Dykman2 . Because $S$ is minimized, the probability of the corresponding trajectory is maximized– a consequence of the WKB approximation. Therefore, all that is needed to find the most probable path to extinction (OP), and $\rho(\mathbb{x},t)$ (Eq.(2) and Eq.(7)) is the appropriate solution of Eqs.(5-6). Such solutions can be determined from boundary conditions, and computed as detailed in Sec.III.

From inspection of Fig.1 we notice that $\rho(\mathbb{x},t)$ is a maximum at the endemic equilibrium ($\mathbb{x}\!=\!\mathbb{x}^{*}$), implying $\partial S/\partial\mathbb{x}\!=\!\mathbb{0}$ and a boundary condition: $\dot{\mathbb{x}}\!=\!\mathbb{0}$, $\dot{\mathbb{p}}\!=\!\mathbb{0}$, $\mathbb{x}\!=\!\mathbb{x}^{*}$, and $\mathbb{p}\!=\!\mathbb{0}$. Second, at the extinct state $(\mathbb{x}\!=\!\mathbb{0})$ the distribution has negative slope, $\partial S/\partial\mathbb{x}\!<\!\mathbb{0}$. Furthermore, $\rho(\mathbb{x},t)$ is approximately time-independent, or quasi-stationary Dykman2 ; Assaf1 ; Assaf2 . In the $WKB$ ansatz, we have $\partial S/\partial t\!=\!H\!=\!0\;\forall t$. Therefore the final boundary condition is $\dot{\mathbb{x}}\!=\!\mathbb{0}$, $\dot{\mathbb{p}}\!=\!\mathbb{0}$, $\mathbb{x}\!=\!\mathbb{0}$, and $\mathbb{p}\!=\!\mathbb{p}^{*}$ Assaf1 ; Schwartz3 ; Dykman . OPs computed with Eqs.(5-6) and the stated boundary conditions are compared with stochastic trajectories ending in extinction for several networks in Fig.(2).

Two important details should be pointed out. Since the distribution is time-independent, and therefore “zero energy”, $S(\mathbb{x})$ is simply the line integral of the momentum along the OP, from Eq.(7). Also, the $\mathbb{p}\equiv\mathbb{0}$ solution of Eqs.(5-6) gives the familiar quenched mean-field equations for SIS model on complex networks. Therefore, the WKB approach generalizes mean-field results to include large fluctuations.

In general, by studying Eqs.(5-6) we can learn how a network’s infection density is coupled to its large fluctuations – together generating the most likely transition sequence through a network leading to extinction. In addition to the distribution of large fluctuations, Eq.(2), an important observable from the above formalism is the geometry of the OP; e.g., specifying the shape of infection density in the different graph positions as a network makes its way from a large epidemic to extinction. Examples are shown in Fig.2 and Fig.5. Another important observable is the average extinction time for a given network, $\left<T\right>$, which is expected to take the form:

[TABLE]

from the assumption that absorption into the extinct state has a rate, or inverse time, proportional to the probability Doering ; Meerson ; Assaf1 ; Schwartz3 . For sufficiently large $S$, the exponential contribution dominates, and therefore $\left<T\right>\!\sim\!e^{S(\mathbb{0})}$ (as demonstrated in Fig.(3) for several networks FN1 ). We note that beyond a theoretical interest, the framework presented can be augmented with control strategies designed to minimize the Action, Eq.(7), thus producing exponential and optimal reductions in the lifetime of epidemics on networks Hindes1 .

II.1 Near threshold behavior

Since the OPs are in a high $2N$-dimensional space, they must be found by numerically solving the two-point boundary value problem in general, given Eqs.(5-6). However, analytic properties can be derived in certain limiting cases, which are useful for guiding intuition and for initializing algorithms (see Sec.(III.1)). An important case discussed in this section is for $\tilde{\beta}$ just above the epidemic threshold, or transcritical bifurcation, $\tilde{\beta}_{c}\!\equiv\!1/\lambda^{(1)}$ – where $\lambda^{(1)}$ is the largest eigenvalue of $A$ Pastor ; Goltsev ; Mata . At this point the endemic and extinct state meet, and below which no long-lived epidemic occurs.

In order to describe the path when $\tilde{\beta}\!\gtrsim\!\tilde{\beta}_{c}$, it is useful to assume that $\tilde{\beta}\!=\!(1+\delta)/\lambda^{(1)}$, with $\delta\ll 1$, and first find the equilibria as functions of $\delta$. Substituting the series $x_{i}^{*}\!=\!\sum_{n}\delta^{n}x_{i,n}$ for each $i$ and $\mathbb{p}\!=\!\mathbb{0}$ into $d\mathbb{x}/d\tau\!=\!\mathbb{0}$, and collecting powers of $\delta$ (e.g., to $\mathcal{O}(\delta^{2})$) gives:

[TABLE]

Furthermore, by decomposing Eqs.(9-10) into the eigenbasis of $A$, $x_{i,n}=\!\sum_{m=1}^{N}G_{n}^{(m)}\eta_{i}^{(m)}$, where $\bm{\eta}^{(m)}$ is the $m$th right eigenvector of $A$ with eigenvalue $\lambda^{(m)}$, and taking the inner product, $\sum_{i}\zeta_{i}x_{i}$, of Eqs.(9-10) with the left eigenvector, $\bm{\zeta}^{(1)}$, we find $\mathbb{x}^{*}$ to $\mathcal{O}(\delta)$. A similar procedures gives $\mathbb{p}^{*}$, when $\mathbb{x}\!=\!\mathbb{0}$:

[TABLE]

(assuming the normalization $\sum_{i}\zeta^{(1)}_{i}\eta^{(1)}_{i}\!=\!1$). Examining Eqs.(11-12), we see that $x_{i}^{*}$ and $p_{i}^{*}$ are proportional to the principal right and left eigenvectors of $A$ near the bifurcation, respectively. In particular, if $\bm{\eta}^{(1)}$ and $\bm{\zeta}^{(1)}$ contain relatively few nodes that are significantly large compared to most others, we expect the infection density and fluctuations to be localized around these nodesGoltsev .

Further insight on the effects of topology near threshold can be gained by considering the Action along the path, Eq.(7). To this end, it is useful to introduce a length parameter, $a\in[0,1]$, so that we can express the coordinates as $x_{i}(a)\!\approx\!x_{i}^{*}(1-a)$ and $p_{i}(a)\!\approx\!p_{i}^{*}a$, where the linear form is the simplest satisfying the boundary conditions to $\mathcal{O}(\delta)$ (see Fig.5 insets). Integrating over the path gives the action near bifurcation, $S(a)\!=\!\sum_{i}\int_{0}^{a}{p_{i}(a^{\prime})da^{\prime}(dx_{i}/da^{\prime})}$:

[TABLE]

We note that Eq.(13) is interesting, since the known expression for the complete graph is generalized by a topological factor that depends on the moments of the centrality distribution. Typically, as the distribution becomes broad, the topological factor in Eq.(13) is reduced, such that the Action differs significantly from the limiting case, $\eta_{i}\!=\!\zeta_{i}\!=\!1/\sqrt{N}\Rightarrow S=N\delta^{2}(a-a^{2}/2)$ Assaf1 . This is intuitive, since for heterogeneous networks, infection is most prevalent around a comparatively small number of nodes, who must recover without reinfection in order for extinction to occur. The effects of heterogeneous eigenvector centrality are explored in more detail in Sec.III.

III SPECIAL SOLUTIONS

In general, the OP is of interest away from threshold. However, since the OP is a heteroclinic connection of Eqs.(5-6), in practice it must be constructed numerically, e.g., through shooting, or quasi-newton methods, etc. For example, the paths shown in Fig.(2) were found from an iterative action minimizing method (IAMM) Lindley2 . In the IAMM, OPs are generated from a least-squares algorithm that minimizes the residuals between Eqs.(5-6) and finite-difference approximations. The boundary conditions specified in Sec.II are used to close the differencing. Often the small $\delta$ limit, Eqs.(11-12), can be used as an initial guess. However the dimension for the minimization is $2Nd$ 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$. Therefore, in practice it is necessary to coarse-grain the network in some way. Two such approaches are discussed in the following sections for networks with large spectral gaps (Sec.III.1) and specified degree distributions (Sec.III.2).

III.1 Large spectral gaps

In general, $A$ can be usefully expanded in terms of its eigenvalues and eigenvectors: $A_{ij}\!=\!\sum_{n=1}^{N}\lambda^{(n)}\eta^{(n)}_{i}\zeta^{(n)}_{j}$. Of particular interest, for strongly connected graphs, $\bm{\eta}^{1}$ and $\bm{\zeta}^{1}$ are positive and unique, and $\lambda^{(1)}$ is equal to the spectral radius, by the Perron-Frobenius theorem. Moreover, in many strongly connected networks of interest, it is the case that $\lambda^{(1)}\!\gg\!\lambda^{(n)}$, with large spectral gaps, and therefore, $A_{ij}\!\approx\!\lambda^{(1)}\eta^{(1)}_{i}\zeta^{(1)}_{j}$. In such cases, $A$ can be coarse-grained along a single dimension as demonstrated below.

Given a large spectral gap, a simple coarse-graining is to bin $\bm{\eta}^{(1)}$, assuming a number of bins, $B$, and a distribution, $f_{b}$, for the number of nodes in a given bin, $b$. In the following, we assume that $A$ is symmetric, so that $\bm{\eta}^{(n)}\!=\!\bm{\zeta}^{(n)}$. In this case, nodes can be ordered according to increasing $\eta_{l}^{(1)}$. The binning procedures follows: starting with the first node, the first bin is filled with nodes sequentially until the number of nodes equals $Nf_{1}$; then, the second bin is filled, etc. Once all nodes are binned, the dimension of Eqs.(5-6) can be reduced by replacing $\eta_{l}^{(1)}$, $x_{l}$, and $p_{l}$ with their bin averages $\forall l\!\in\!b$: $\eta_{b}\!\equiv\!\sum_{l\in b}\eta_{l}^{(1)}\!/(Nf_{b})$, $x_{b}\!\equiv\!\sum_{l\in b}x_{l}/(Nf_{b})$, $p_{b}\!\equiv\!\sum_{l\in b}p_{l}/(Nf_{b})$. This gives the following approximations to Eqs.(5-6) with reduced dimension $2B$:

[TABLE]

A final requirement is needed to ensure that the binned and original system have the same bifurcation point. We choose to renormalize $\eta_{b}$ so that $\sum_{b}\eta_{b}^{2}f_{b}N\!=\!\sum_{i}\eta_{i}^{2}\!=\!1$.

The above procedure was applied to three networks (considered in Figs.2-3), where the bin distribution was assumed to be uniform for simplicity, $f_{b}\!=\!1/B$; $B$ was chosen large enough so that the binned centralities closely matched the original (as in Fig.4), but not so large to preclude using $200\!-\!1000$ discretization points along the $OP$ in the IAMM. The first network in Fig.2(a) is a weighted Barabási-Albert graph (WBA) with $N\!=\!500$ and initial degree for each node, $m\!=\!7$ Albert . Every link was given a random weight, independently drawn from a uniform distribution over the range $[0,10]$ after the network was generated from the standard Barabási-Albert algorithm ($\lambda^{(1)}\!=\!122,\lambda^{(2)}\!=\!58.6$). The second network in Fig.2(b) is a high-resolution American high school contact network (HS) with $788$ individuals and links representing close proximity interactions during the course of a day (measured using wireless motes Salathe ). Weights associated with links correspond to contact durations ($\lambda^{(1)}\!=\!6715,\lambda^{(2)}\!=\!4882$). Binning results for the WBA and HS networks are shown in Fig.4 with $B\!=\!55$ and $B\!=\!30$, respectively. The third network was generated from a configuration model, $CM$, with $N\!=\!1000$ and a degree distribution, $g(k)\!=\!k^{-2.5}\!/\!\sum_{k^{\prime}=20}^{500}k^{\prime-2.5}$, where $k$ is the number of links for a given node ($\lambda^{(1)}\!=\!80.9,\lambda^{(2)}\!=\!16.3$)Newman2 ; $B\!=\!55$ for the CM network’s OP computation. A related method for computing OPs through networks with specified degree distributions is discussed in Sec.III.2, and was applied to the networks in Fig.2(c)-(d).

III.1.1 Scaling with broad centrality distribution

As $\tilde{\beta}$ is increased above $\tilde{\beta}_{c}$, infection increases across the network. In particular, infection density can become high even at low centrality nodes. In this case the path to extinction through a network is more complex as its global dynamical structure becomes apparent. Example optimal paths for a WBA network when $\tilde{\beta}\!\gg\!\tilde{\beta}_{c}$ are shown in Fig.5. We can see that a multi-step structure is visible in the relative change of infection density at different graph positions, which can be compared with the insets ($\tilde{\beta}\!\gtrsim\!1/\lambda^{(1)}$) Hindes1 . Though the path has a more complicated form away from threshold, some characteristic scalings can be captured in this region of parameter space for networks with heterogeneous eigenvector centralities (e.g., $f_{b}\sim\eta_{b}^{-\gamma}$) and where the large spectral gap approximation holds. Our approach in this section is to study the unstable and stable linear modes of $(\mathbb{x},\mathbb{p})$ near the endemic and extinct states, respectively, for such networks. These modes approximate the OP (a heteroclinic connection) near the equilibria, and are useful for describing how large fluctuations depend on centralityFN6 .

With this end in mind, we consider the dynamics of $(x_{i},p_{i})=(x_{i}^{*}+\epsilon_{i}^{o},\mu_{i}^{o})$ and $(x_{i},p_{i})=(\epsilon_{i}^{in},p_{i}^{*}+\mu_{i}^{in})$, for small $\bm{\epsilon}$ and $\bm{\mu}$, given $A_{ij}\approx\lambda^{(1)}\eta^{(1)}_{i}\eta^{(1)}_{j}$ (below, we drop the superscript (1) in $\bm{\eta}$ and $\lambda$ for convenience). Similar to Sec.II.1, when $\tilde{\beta}\!\gtrsim\!\tilde{\beta}_{c}$, it is straightforward to show that $\epsilon_{i}^{o},\epsilon_{i}^{in},\mu_{i}^{o}$ and $\mu_{i}^{in}$ are simply proportional to $\eta_{i}$. Fig.6 shows centrality scalings for the principal linear eigen-modes of Eqs.(5-6) near the equilibria. The upper dashed lines demonstrate the predicted scaling $\epsilon_{i}^{in}/\epsilon_{j}^{in}\!\sim\!\eta_{i}/\eta_{j}\!\sim\!\mu_{i}^{o}/\mu_{j}^{o}$ for a WBA network where the dark blue/red curves correspond to $\tilde{\beta}$ increasingly close to threshold ($\epsilon_{i}^{o}$ and $\mu_{i}^{in}$ scale similarly in this region). However, as $\tilde{\beta}$ is increased, shown in light blue/red, we can see that the scaling changes significantly Hindes1 .

In order to understand the change in scaling as $\tilde{\beta}$ is increased, we first consider the equilibria $x_{i}^{*}$ and $p_{i}^{*}$. Given the large spectral gap assumption, we find a simple form for each that is dependent on two parameters, $X$ and $P$:

[TABLE]

satisfying: $X\!=\!\sum_{j}\eta_{j}x^{*}_{j}$ and $P\!=\!\sum_{j}\eta_{j}e^{p_{j}^{*}}$, where $\left<\eta\right>$ is the average eigenvector centrality Ott2 . In particular, by assuming that infection densities are high in the endemic state at most graph positions, i.e., $\tilde{\beta}\lambda\eta_{i}N\!\left<\eta\right>\!\gg\!1$, then $X\!\approx\!N\!\!\left<\eta\right>\!-\!1/[\tilde{\beta}\lambda\!\left<\eta\right>]$ and $P\!\approx\!1/[\tilde{\beta}\lambda\!\left<\eta\right>].$ Substituting these approximations into the linearized Eqs.(5-6) allows us to determine the dependence of the eigen-modes, $(\epsilon_{i}^{o}(t),\mu_{i}^{o}(t))\!=e^{t\sigma^{o}}\!(\epsilon_{i}^{o},\mu_{i}^{o})$ and $(\epsilon_{i}^{in}(t),\mu_{i}^{in}(t))\!=e^{t\sigma^{in}}\!(\epsilon_{i}^{in},\mu_{i}^{in})$, on $\eta_{i}$ near the equilibria.

Since infection densities are high near the endemic state away from threshold, we expect the most well connected nodes to be quickly reinfected after recovery, as compared to nodes that are less well connected. Therefore, we expect the OP out of the endemic state to correspond with an initial decrease in infection at low centrality positions. The scaling for this initial step is determined by the eigensolution of the linearized Eq.(6) at $(x_{i}^{*},0)$:

[TABLE]

In particular, $\sigma^{o}$ is positive and grows from zero with $\tilde{\beta}\!>\!\tilde{\beta}_{c}$ FN2 . When $\tilde{\beta}\lambda\eta_{i}N\!\left<\eta\right>\!\gg\!1$ and $f_{b}\!\sim\!\eta_{b}^{-\gamma}$ for large $\eta$, the summation in Eq.(18) $\sim\!\int\!\eta^{-\gamma+1}d\eta/(\eta-\text{const.})$, and converges in the limit of large maximum centrality, $\eta_{max}$, when $\gamma\!>\!2$. This implies that $\sigma^{o}$ does not depend sensitively on $\eta_{max}$ in this region and therefore we can consider the limit of large centrality in Eq.(19). By inspecting $\mu_{i}^{o}$ for large $\eta_{i},$ we see that it tends to a constant, i.e., $\mu_{i}^{o}/\mu_{j}^{o}\sim 1$, (since the sum over $j$ is $i$-independent) in good agreement with numerical solutions of Eqs.(18-19) away from threshold– shown in Fig.6(b)(light red). Interestingly, $\mu^{o}_{i}$ becomes largest for small $\eta_{i}$ (light red) and increases quickly to a constant for large $\eta_{i}$. Since the momenta are nearly equal across nodes in the network, the Action’s derivatives w.r.t. infection density are nearly equal across nodes. A similar procedure gives the scaling $\epsilon_{i}^{o}/\epsilon_{j}^{o}\sim\eta_{j}/\eta_{i}$, which is found by expanding the linearized Eq.(5) in $1/\tilde{\beta}\lambda\eta_{i}N\!\left<\eta\right>$, e.g., $x_{i}^{*}\!\approx\!1-1/\tilde{\beta}\lambda\eta_{i}N\!\left<\eta\right>$ (Sec.VI). Therefore, the infection density at a given node decreases inversely proportional to its eigenvector centrality, i.e, the reciprocal scaling of the OP near threshold.

Following the same approach, the scaling of the OP near the extinct state is determined by the eigen-solution of the linearized Eq.(5) at $(0,p_{i}^{*})$:

[TABLE]

In particular, $\sigma^{in}$ is negative, decreases from zero with $\tilde{\beta}\!>\!\tilde{\beta}_{c}$, and is similarly insensitive to large $\eta_{max}$. Hence, taking the limit of large $\eta_{i}$, given e^{-p_{i}^{*}}\!\approx\!1\!+\!\beta\lambda\eta_{i}N\!\left<\eta\right>\![1\!-\!1\big{/}\beta\lambda\eta_{i}N\!\left<\eta\right>], implies $\epsilon_{i}^{in}/\epsilon_{j}^{in}\sim\eta_{j}/\eta_{i}$ in Eq.(21)– as found near the endemic state and shown in Fig.6(b)(light blue). However, in contrast to the behavior near $x_{i}^{*}$, we find that $\mu_{i}^{in}$ increases with $\eta_{i}$, for small $\eta_{i}$, before reaching a constant, $\mu_{i}^{in}/\mu_{j}^{in}\sim 1$ (see details in Sec.VI and Fig.9). The scalings near the extinct state imply that the last segment of the OP is coincident with a final recovery of residual infections at low-centrality nodes, while the momentum is largest at high centralities.

Putting the scalings near the equilibria together, we can infer that infection density decreases rapidly in high-centrality graph positions at a boundary layer between the endemic and extinct states – since we have shown that their change is small compared to low centralities near the equilibria Hindes1 ; FN3 . This can be seen in Fig.5(a), where projections of the OP into high and low-centralities, on the $x$ and $y$ axes respectively, show a characteristic pattern in which segments with large horizontal slope occur between two segments with large vertical slope.

III.2 Degree distributions

In addition to understanding the OP for a given network defined by $A$, it is useful to understand the qualitative structure of paths and Actions for networks with similar statistical properties Pastor ; Hindes1 ; Dorogovtsev . A popular approach is to consider networks with a specified distribution for the fraction of nodes with $k$ links, $g(k)$ (where $k$ is called the degree), which is the focus of this section. Often, additional information is stipulated, such as a degree-correlation function – typically in the form of a specified probability that a link starting from a node with degree $k$ leads to a node with degree $k^{\prime}$, $o(k^{\prime}|k)$ Vespignani ; Pastor2 ; Pastor3 . OPs for networks with such properties can be found by approximating $A$ given these distributions, and substituting into Eqs.(5-6).

As is customary, we replace $A_{ij}$ by its expectation value in the ensemble of simple networks with $g(k)$ and $o(k^{\prime}|k)$, which is called the annealed network approximation. In particular, $A_{ij}$ is approximated by the probability that nodes $i$ and $j$ are connected, A_{ij}\!\approx\!o(k_{j}|k_{i})k_{i}\big{/}Ng(k_{j}), or the probability that node $i$ is connected to any node with degree $k_{j}$ along a single link, multiplied by the number of possible links, and divided by the number of nodes with degree $k_{j}$ Dorogovtsev ; Pastor ; Hindes1 . Note that for link consistency, $A_{ij}=A_{ji}$, the distributions must satisfy the constraint: $ko(k^{\prime}|k)p(k)=k^{\prime}o(k|k^{\prime})p(k^{\prime})$ Pastor3 . With this substitution for $A_{ij}$ into Eqs.(5-6), Hamilton’s equations depend on the density of infection for nodes with the same degree $k$, $x_{k}$, and their momentum, $p_{k}$:

[TABLE]

Notably, Eqs.(22-III.2) reduce to the heterogeneous mean-field dynamics for networks when $p_{k}\!\equiv\!0$; $p_{k}\!\neq\!0$ entails extinction in degree-correlated topologies with dimension of $(\mathbb{x},\mathbb{p})$ equal to twice the number of degree classes. The analysis and results for degree-distributed networks are analogous to Sec.II.1-Sec.III.1.1. For example, for degree-distributed networks the familiar proportionality of the Action on the number of nodes in $A$, is found from Eq.(7) Dykman2 :

[TABLE]

Moreover, with the appropriate substitution of the largest eigenvalue, $\lambda$, and the corresponding right eigenvector, $v_{k}$, of $ko(k^{\prime}|k)$ in Eq.(13) Pastor3 , we find the Action at the extinct state for degree-correlated networks near the epidemic threshold:

[TABLE]

where $\left<v^{n}\right>\!=\!\sum_{k}g(k)v_{k}^{n}$.

Extinction paths and times for two example networks are shown in Figs.2-3 Hindes1 ; Hindes4 . The networks have a bimodal degree-distribution with positive $(PC)$ and negative ($NC$) degree-correlations (Figs.2(c) and (d) respectively), where positive implies an increased probability relative to an uncorrelated network for nodes with similar degree to share an edge. Correlated bimodal networks can be constructed in a straightforward manner as detailed in Sec.VII. Fig.2 (c)-(d) shows OPs for example parameters computed from Eqs.(22-III.2), and projected into the densities of infected low and high-degree nodes. Qualitatively, we can see that OP projections into infection densities are significantly closer to lines with unit slope (which is the case for uncorrelated networks with small variance in $k$) in the $NC$ case (d), than for the $PC$, (c).

The change in the OP’s shape with correlations suggests a reduction/enhancement of the effects of network heterogeneity with negative/positive correlations. For positive correlation, infection is more prevalent around high-degree nodes. This is reflected in the principal eigenvectors of $ko(k^{\prime}|k)$ for the two examples, where the low-degree component is $5.4$ times greater in the $NC$ (parameters in Fig.2). In fact the topological factor, $\left<v^{2}\right>^{3}\!\!/\!\left<v^{3}\right>^{2}$, is $2.2$ times greater for the $NC$, given the same $N$ and distance to bifurcation, $\delta\!=\!\tilde{\beta}\lambda\!-\!1\!\gtrsim\!0$. Therefore we expect the probabilities for large fluctuations to be smaller by the same power and extinction times to be larger by the same power. Equivalently, if comparing fixed extinction time, the $PC$ must be taken to larger $\tilde{\beta}\lambda$ and/or $N$. This is demonstrated in Fig.3, where the largest times shown correspond to $\tilde{\beta}\lambda\!=\!1.9$ and $N\!=\!400$ for the $NC$, and $\tilde{\beta}\lambda\!=\!2.8$ and $N\!=\!500$ for the $PC$.

The above example raises an interesting question of how fluctuations and extinction times vary with statistical properties in a network, such as degree-heterogeneity, which can be anticipated from the network Action. A more realistic class of heterogenous networks have power-law degree distributions, $g(k)\!\sim\!k^{-\gamma}$, where the level of degree-heterogeneity grows with decreasing $\gamma$. Fig.7 shows the predicted Actions at the extinct state as a function of $\gamma$ for truncated and uncorrelated, $o(k|k^{\prime})\!=\!kg(k)/\!\left<k\right>$, power-law distributions with several fixed distances to threshold. Interestingly, for such networks we can see that Actions vary as much as $60\%$ when ($\tilde{\beta}\lambda\!=\!2$), with broader distributions resulting in significantly smaller Actions, and therefore exponentially larger probabilities of large fluctuations and exponentially smaller extinction times.

The Action curves were found by solving Eqs.(22-III.2) with the boundary conditions specified in Sec.II, and computing Eq.(24). In addition to computation, the lower black-curve in Fig.7 gives the analytic scaling near threshold, found from Eq.(25), by substituting the eigenvector for uncorrelated random networks, $v_{k}\!=\!k\!\Rightarrow S\!=\!N\delta^{2}\!\left<k^{2}\right>^{3}\!\!/\!\left<k^{3}\right>^{2}$ Pastor3 ; Hindes1 . For the computed curves, it was useful to reduce the dimension for the IAMM by binning the distribution $o(k|k^{\prime})\!=\!kg(k)/\!\left<k\right>$ with a similar procedure as Sec.III.1. Our approach, was to select a small bin width for $kg(k)/\!\left<k\right>$ (e.g., 0.015), and sequentially add degree classes to a bin, starting with the smallest $k$ and first bin, until the sum of $kg(k)/\!\left<k\right>$ over $k$ in a given bin equaled or exceeded the bin width. Then, the next bin was filled with the same bin width, etc. In the final step, degrees in Eqs.(22-III.2) were replaced by their bin’s average and $o(k^{\prime}|k)$ by the sum over $kg(k)/\!\left<k\right>$ in each bin Hindes1 .

III.2.1 System-size scaling for modest N

Another interesting feature of extinction times in power-law networks concerns their scaling with system-size when there is no truncation in $g(k)$. It is known for degree-homogeneous networks, such as simple-complete or Erdős-Rényi graphs, that the Action scales linearly with the system size Assaf1 ; Lindley . Below we show that near threshold for modest $N$, a range of scalings are possible depending on the exponent, $\gamma$. As indicated above, the Action near threshold depends on a topological factor that is a function of the moments of $g(k)$, which can depend on $N$.

Here we continue to use the annealed network approximation, though for very large $N$ this is known to break down for random networks with unbounded degree as localization effects become important Goltsev . Therefore, we restrict ourselves to $N$ and minimum degree, $k_{min}$, such that $\lambda\approx\left<k^{2}\right>\!/\!\left<k\right>$. When $\gamma\!>\!4$, $\left<k^{3}\right>$ is finite for power-law networks, i.e., independent of $N$ for large $N$, and thus $S\!\sim\!N$ near threshold, $\tilde{\beta}\!\gtrsim\!\left<k\right>\!/\!\left<k^{2}\right>$ Vespignani ; Pastor2 – though higher-order terms may be $N$-dependent for $\delta\!\ll\!1$, we expect them to grow more slowly than $N$ FN5 . On the other hand, if $\gamma\!<\!4$, $\left<k^{3}\right>$ is a function of the maximum degree, $k_{max}$, which follows a simple scaling: $k_{max}\sim k_{min}N^{1/[\gamma-1]}$, for a finite network with minimum degree $k_{min}$ Cohen . The customary approach is to approximate the statistical moments of $g(k)$ given $k_{max}$, allowing one to find the scaling of $S$ with $N$. For example, when $k_{min}\!\gg\!1$, the discrete sum, $\left<k^{3}\right>\!=\!\sum_{k}g(k)k^{3}\approx\tilde{C}\int_{k_{min}}^{k_{max}}k^{3-\gamma}dk$. Introducing $\gamma\!=\!3\!+\!\alpha$ with $\alpha\!\in\!(0,1)$, we get $\left<k^{3}\right>\!\approx\!\tilde{C}k_{min}^{1-\alpha}N^{[1-\alpha]/[2+\alpha]}/[1-\alpha]$, where $\tilde{C}\!=\!k_{min}^{2+\alpha}[2+\alpha]$ (from normalization of $g(k)$). Computing the moments of $g(k)$ in this way, gives the Action at the extinct state to $\mathcal{O}(\delta^{2})$ from Eq.(24):

[TABLE]

The above suggests that in the heterogeneous mean-field approximation, the $\mathcal{O}(\delta^{2})$ contribution to the Action can increase sub-linearly in $N$ for $\gamma\!\in\!(3,4)$ near threshold Assaf3 as suggested in Fig.8. However for very large networks, and no truncation in $k_{max}$, eventually $\lambda\!\sim\!\text{max}\{\sqrt{k_{max}},\left<k^{2}\right>\!/\!\left<k\right>\}\!\gg\!1$, and the analysis presented is no longer valid, including the expansion in $\delta$. Moreover, there is some evidence for multiple epidemic thresholds in networks with unbounded $k_{max}$ as $N\rightarrow\infty$ Mata2 . Since such issues are not yet resolved, we leave the description of extinction in very large networks with unbounded degree distributions, and the crossover between localized and delocalized extinction for future study.

IV CONCLUSION

This work dealt with the extinction of long-lived endemic states above epidemic thresholds on static finite networks with infection dynamics given by the stochastic SIS model. The optimal path to extinction (OP), the distribution of large fluctuations, and the average extinction time were computed by combining mean-field and WKB-approximation techniques. The path-based formalism presented enabled us to predict extinction in general networked populations, and extract several of its intriguing signatures in complex topologies, including the multistep scaling of the OP in networks with heterogeneous eigenvector centrality, as well as an increase in the probability of large fluctuations with increased topological heterogeneity. Although theoretical in nature, the generality of our approach allowed us to consider several applications, including weighted empirical and degree-correlated topologies.

Though the results show good qualitative and quantitative agreement with Monte-Carlo simulations in both real and synthetic networks, improved accuracy can be achieved in a straightforward manner by following our synthesized prescription, namely: Using as an ansatz in a network’s master equation the exponential function of an Action (typically requiring some accurate mean-field approximation), and taking a large system-size limit. The result is a Hamilton-Jacobi equation that generates a dynamical system with twice the dimension of the mean-field. The OP can be found by solving the two-point boundary value problem of Hamilton’s equations of motion beginning at an endemic state and ending at an extinct state, which define OP endpoints. Thus the theory changes the stochastic analysis of large fluctuations in networks to one that may be analyzed using a deterministic formalism whose zero-fluctuation limit is a mean-field theory. Furthermore, our approach can be more generally applied to other questions concerning noise and network dynamics, such as epidemic extinction in adaptive networks, switching in social networks, network inference in the presence of large fluctuations, and optimal control of networks with fluctuating dynamics Hindes1 ; Motter .

V ACKNOWLEDGMENTS

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 C. R. Myers and M. Assaf for useful discussions.

VI APPENDIX A

As described in Sec.III.1.1, $\epsilon_{i}^{o}$ and $\mu_{i}^{in}$ satisfy the linearized Eqs.(5-6). When $\tilde{\beta}\lambda\eta_{i}N\!\left<\eta\right>\!\gg\!1$, the approximate linear systems are:

[TABLE]

[TABLE]

Since the sums in Eqs.(27-28) are independent of $i$, the limit of large $\tilde{\beta}\lambda\eta_{i}N\!\left<\eta\right>$, gives: $\epsilon_{i}^{o}/\epsilon_{j}^{o}\sim\eta_{j}/\eta_{i}$ and $\mu_{i}^{in}/\mu_{j}^{in}\sim 1$. The latter can be seen in Fig.(9). However, as $\tilde{\beta}\!\gg\!\tilde{\beta}_{c}$ the continuous spectra, $\sigma^{o}_{i}$ and $\sigma^{in}_{i}$, for large $N$ of the linearized Eqs.(5-6) become relevant: $\epsilon_{i}^{in},\mu_{i}^{o}\sim\delta{(\eta-\eta_{i})}$, and

[TABLE]

This occurs as the denominators of Eq.(18) and Eq.(20) approach zero, and the single-mode analysis of Sec.III.1.1 is invalid. As a consequence, for very large $\tilde{\beta}$, the relevant modes directing the OP to extinction near the equillibria are extremely localized around low-centrality nodes.

VII APPENDIX B

Correlated bimodal networks can be constructed as follows. We assume that a fraction, $p$, of the network has high-degree near $k_{2}$ while the remaining nodes have low-degree near $k_{1}$. To build such networks, high-degree nodes are connected to each other with probability $k_{2}^{2}/[N\!\left<k\right>]\!+w$, where $w$ measures the assortativity above the uncorrelated construction and $\left<k\right>\!=\!(k_{1}(1-p)\!+\!k_{2}p)$. On the other hand, high and low-degree nodes are connected with probability $k_{1}k_{2}/[N\!\left<k\right>]-w$, and low-degree nodes are connected with probability $k_{1}^{2}/[N\!\left<k\right>]\!+\!w^{\prime}$ – where $w^{\prime}$ is determined from the link-consistency constraint. In this way the degree distribution has two peaks centered around $k_{1}$ and $k_{2}$ as $N\rightarrow\infty$, and Eqs.(22-III.2) can be used to capture the OP and average extinction times assuming two degree classes with: $o(k_{2}|k_{2})\!=\!w+k_{2}p/\!\left<k\right>$, $o(k_{1}|k_{2})\!=\!-w+k_{1}(1-p)/\!\left<k\right>$, $o(k_{1}|k_{1})\!=\!-w^{\prime}+k_{1}(1-p)/\!\left<k\right>$, and $o(k_{2}|k_{1})\!=\!-w^{\prime}\!+\!k_{2}p/\!\left<k\right>$.