Physical Review Letters
American Physical Society
Epidemic Extinction and Control in Heterogeneous Networks
DOI 10.1103/PhysRevLett.117.028302, Volume: 117, Issue: 2,




We consider epidemic extinction in finite networks with a broad variation in local connectivity. Generalizing the theory of large fluctuations to random networks with a given degree distribution, we are able to predict the most probable, or optimal, paths to extinction in various configurations, including truncated power laws. We find that paths for heterogeneous networks follow a limiting form in which infection first decreases in low-degree nodes, which triggers a rapid extinction in high-degree nodes, and finishes with a residual low-degree extinction. The usefulness of our approach is further demonstrated through optimal control strategies that leverage the dependence of finite-size fluctuations on network topology. Interestingly, we find that the optimal control is a mix of treating both high- and low-degree nodes based on theoretical predictions, in contrast to methods that ignore dynamical fluctuations.

Hindes and Schwartz: Epidemic Extinction and Control in Heterogeneous Networks

The extinction of epidemics in finite networks is an important topic in population dynamics [1,2]. Though many factors may contribute, such as environmental changes and social behavior, it has been demonstrated, and rigorously proven for finite populations, that internal fluctuations in a system’s dynamics can organize in such a way to induce a large fluctuation along a most probable, or optimal, path to extinction [3–5]. Such fluctuations to infection-free states have been studied extensively inwell-mixed systems, including the role of vaccination and treatment programs in reducing the average time to extinction [3,6]. Similarly, the most probable extinction paths have been found in networks with homogeneous degree, but the behavior appears to be independent of network topology, as in the well-mixed limit [7].

Somewhat separately, much work has been done in characterizing the deterministic dynamics, epidemic threshold, outbreak size distributions, small fluctuations, localization, and phases of epidemics in complex networks [8–14]. Only very recently has there been progress in understanding the interplay between stochastic noise and network dynamics that can lead to large fluctuations and switching between states [15]. However, very little is known about how internal noise inherent to epidemic extinction pertains in heterogeneous networks having vastly differing topologies.

In this Letter we construct and analyze the most probable path through heterogeneous networks to extinction. Novel in this work is that we show the path has two primary forms, close to and far from the epidemic threshold. In the latter, we demonstrate an interesting multistep structure in which low-degree node infections decrease first, followed by a quick, nearly complete extinction in high-degree nodes, and finishing with a low-degree extinction. The approach is then used to design a novel targeted optimal treatment strategy that can exponentially reduce extinction times, but does not trivially treat the most well connected nodes, minimize the epidemic size, or maximize the number of treatments. Instead, the optimal control minimizes the “action” associated with a transition to extinction with respect to the network topology. Such controls that manipulate finite-size fluctuations inherent in contact processes based on topology are also novel in the study of complex networks [6,15].

To understand how extinctions depend on topological heterogeneity, we consider the stochastic susceptible-infectious-susceptible (SIS) model on uncorrelated random networks with a given degree distribution gk, where the degree k is the number of links of a node. Simple graphs with N nodes can be generated from gk in several ways, for example, with a configuration model network (CMN) [16]. Such networks are usefully represented by an adjacency matrix A, where Aij is 1 if nodes i and j are linked, and 0 otherwise. In this representation a network’s SIS dynamics is captured by the states and transitions of its nodes; e.g., node i is either infected, denoted νi=1, or susceptible, νi=0, and changes its state νi:01 with probability per unit time β(1-νi)jAijνj, and νi:10 with probability per unit time ανi, where β and α are known as the infection and recovery rates, respectively.

In order to describe the dynamics given these reactions, it is useful to approximate Aij with its expectation value in an ensemble of networks Aij, which for uncorrelated networks takes the form Aij=kikj/(Nk) in the limit of large N. This is known as the “annealed” network approximation and represents a mean field for heterogeneous networks [8–10,14], though other techniques give similar results to those shown here [17]. Given this form, the state of the network can be described by the number of infected nodes with degree k, Ik, which has corresponding reactions and rates: IkIk+1 with rate of infection wkinf(I)=βk(Nk-Ik)kkIk/(Nk), and IkIk-1 with recovery rate wkrec(I)=αIk, where Ik=Ik1,Ik2,,Ikmax, and Nk=gkN.

Since the SIS model is stochastic, associated with all network states is a probability distribution ρ(I,t), which for networks having a general degree distribution satisfies an approximate master equation

where 1k=0k1,0k2,,1k,,0kmax. We are interested in the behavior of Eq. (1) for large, but finite, networks. As customary, we assume N is large and take the leading order in a 1/N expansion [3,18]. This is similar to the Wentzel-Kramers-Brillouin (WKB) ansatz of quantum mechanics, where 1/N plays the role of Planck’s constant in Schrödinger’s equation. In accordance with WKB, by writing ρ(I,t)=e-NS(x,t), where x=I/N, and taking the leading order in 1/N, or wk(I±1k)wk(I) and ρ(I±1k,t)e-NS(x)eS/xk, we find a Hamilton-Jacobi equation, (S/t)+H[x,(S/x)]=0, where S and H are called the action and Hamiltonian, respectively. As in classical mechanics, the Hamiltonian is a function of the coordinate x and its conjugate momentum p=S/x:
In this context, momenta behave as fluctuations on x—describing both size and direction.

It is convenient to analyze Eq. (2) by solving the canonical equations of motion, x˙k=H/pk, p˙k=-H/xk, in terms of the fraction of each degree class infected yk=xk/gk, the ratio β/α=β˜, and the rescaled time τ=αt:


Of interest are particular solutions of Eq. (3) which correspond to network trajectories that remain near an endemic state for some time, and then decay into an extinct state, with no more infectious nodes. When the two states are well separated, the distribution ρ(I,t) is quasistationary, or in the WKB ansatz, (S/t)=H=0:

This suggests that we look for solutions of Eq. (3) in the form of heteroclinic paths connecting two saddle-point equilibria: from an endemic fixed point, yk*=1/[1+1/(Yβ˜k)], pk=0, to extinction with nonzero momentum, yk=0, pk*=-ln[1+β˜k(1-P)][3,5,19]. The functions Y and P depend on β˜, with (Y,P)(0,1) as β˜k2/kR01, and (Y,P)(1,0) as β˜. Importantly, because such paths extremize their action, they extremize their probability, and therefore correspond to most probable paths through a network to extinction [18,20]. Figure 1 shows comparisons between projections of prehistory trajectories to extinction from stochastic simulations and optimal paths of Eq. (3) for several network configurations computed with the iterative action minimizing method [21].

In general, the average time to reach extinction T will depend on β˜ and network properties in complicated ways [22,23]. However, for sufficiently large N, the transition is an exponential process with a rate proportional to the probability, and therefore TeNS. The network action (4) is thus central to understanding the dependencies of extinction times on network topology.

Qualitatively, we find two important parameter regions to consider. First, close to the epidemic threshold when R0-10 (which we call “weak”), the paths to extinction are approximately linear from (yk*,0) to (0,pk*). The explicit form can be seen by expanding the equilibria in powers of R0-1, which gives to first order yk*kk2(R0-1)/k3=-pk*, implying that the endemic state and momentum at extinction are simply proportional to degree when the infection is weak. This is intuitive since in the weak limit high-degree nodes drive the epidemic, when only their local reproductive numbers are sufficient to spread infection, (β˜k>1), and therefore must recover disproportionately without reinfection in order for extinction to occur. Paths near the weak limit can be seen in Figs. 1(b) and 1(d) where R0=1.6 and R0=2.0, respectively, and in Fig. 2(a) (red).

Density (unnormalized) of 1000 simulations projected into the fraction of infected high-degree (H) and low-degree (L) nodes. Predicted paths are shown in blue from the endemic state (*) to extinction (∘). (a) A network with Aij=kikj/(N⟨k⟩), N=300, and β˜=0.096, and with two degree classes, ki∈{5,50}; k=50 nodes occupying 10% of the network. (b) A corresponding CMN. (c) A CMN with N=350, β˜=0.092, and gk=e-1616k/k! where high-degree nodes have 16≤k≤18, and low-degree have 13≤k≤15. (d) A CMN with N=600, β˜=0.038, and gk=k-2.5/∑k′=10300k′-2.5, where high-degree nodes have 70≤k≤235, and low-degree have 10≤k≤12.
FIG. 1.
Density (unnormalized) of 1000 simulations projected into the fraction of infected high-degree (H) and low-degree (L) nodes. Predicted paths are shown in blue from the endemic state (*) to extinction (∘). (a) A network with Aij=kikj/(Nk), N=300, and β˜=0.096, and with two degree classes, ki{5,50}; k=50 nodes occupying 10% of the network. (b) A corresponding CMN. (c) A CMN with N=350, β˜=0.092, and gk=e-1616k/k! where high-degree nodes have 16k18, and low-degree have 13k15. (d) A CMN with N=600, β˜=0.038, and gk=k-2.5/k=10300k-2.5, where high-degree nodes have 70k235, and low-degree have 10k12.
(a) Projections of the optimal paths for the truncated power law [see Fig. 1(d)] shown for increasing R0∈[1.1,5.1] (red→black, weak→strong) in steps of 0.5, compared with the path into the endemic state for R0=5.1 (green). Arrows indicate direction in time. (b) Projections into yk for the same distribution with R0=9 and 17 bins [24], shown for bins with increasing k{k=12,14≤k≤15,18≤k≤20,24≤k≤27,34≤k≤41,53≤k≤69,97≤k≤143,236≤k≤300} (blue→black), and compared with the predicted scaling for the highest bin (dashed lines).
FIG. 2.
(a) Projections of the optimal paths for the truncated power law [see Fig. 1(d)] shown for increasing R0[1.1,5.1] (redblack, weakstrong) in steps of 0.5, compared with the path into the endemic state for R0=5.1 (green). Arrows indicate direction in time. (b) Projections into yk for the same distribution with R0=9 and 17 bins [24], shown for bins with increasing k{k=12,14k15,18k20,24k27,34k41,53k69,97k143,236k300} (blueblack), and compared with the predicted scaling for the highest bin (dashed lines).

Also for weak infection, the action along the path from Eq. (4) is therefore

which depends on the distance from the epidemic threshold and a nontrivial topological factor that generally decreases with increased broadness in the degree distribution. In contrast, in the well-mixed limit (corresponding to the simple complete graph) the action only depends on R0[3,22]; the predicted reduction in extinction times with topological fluctuations is intuitive, since for very heterogeneous networks only a small fraction of highly connected nodes must recover without reinfection, compared with most nodes in networks where nodes are topologically similar.

On the other hand if most nodes can propagate infection, β˜k1 (which we call “strong”), then the interplay between degree classes and the path to extinction is more complicated as the global dynamical structure of the path becomes apparent. However, we find that a limiting form emerges when comparing the dynamics of low- and high-degree nodes by which the path can be described in multiple steps.

Since in the strong limit most nodes will be infected in the endemic state, it is very improbable that high-degree nodes can recover without being reinfected, and thus we expect infection must first disproportionately decrease in low-degree nodes. We can extract the form of this step, by analyzing the unstable eigenmode of the endemic equilibrium (yk,pk)=(yk*+εk(1),μk(1)) to linear order and studying the asymptotic scaling of (εk(1),μk(1)) for large β˜k. Inserting these assumptions into Eq. (3), we find that (εk(1),μk(1)) must satisfy an eigenvalue equation for the rate λ(1):


Since λ(1) and the sum are k independent, it must be that for large β˜k, we have εk(1)/εk(1)k/k [the relative decrease in infection shown in Fig. 2(b) (1)], with relative momenta initially tending to a constant.

In the second step [Fig. 2(b) (2)], the small buildup of momenta for high-degree nodes becomes rapid as the k-dependent contribution to p˙k approaches a maximum along the path. In analogy with mechanics, this can be thought of as the network’s contribution to the “force” on yk, which near its maximum quickly “pushes” -pk from near zero to its maximum -pk*. On the other hand, since yk and pk decreases along the path, by inspecting y˙k, we can find an upper bound for -y˙k, max[-y˙k]<max[-yke-pk]<-yk*e-pk*, because the k-dependent contribution to y˙k is positive. Since ykyk* until -pk differs from zero, as the force approaches its maximum, -y˙k can be approximated by the upper bound, giving the scaling for large β˜k: εk(2)/εk(2)k/k.

In the last step [Fig. 2(b) (3)], we expect to have a final decrease in low-degree node infections in a background of very small numbers of infected high-degree nodes, since the latter were rapidly depleted in the second step. The scaling can be found by analyzing the extinct state’s stable eigenmode (yk,pk)=(εk(3),pk*+μk(3)), which gives an eigenvalue equation for the rate λ(3)

which in the limit of large β˜k implies εk(3)/εk(3)k/k. Examples are shown in Figs. 2(a) (black) and 2(b) (black); in the former, the strong limit scaling starts to be visible for R03. Also, Fig. 2(a) shows the significant qualitative difference between the epidemic path toward the endemic state (green), which has been well studied, and the optimal path to extinction [8].

In addition to a theoretical interest in the geometry of the optimal path through a network, it is also practically interesting, because extinction times scale exponentially with the action [3]. Since the action depends nontrivially on network topology, e.g., Eq. (5), we suggest exploiting topology as a basis for optimal epidemic control strategies in finite networks, with the goal of minimizing extinction times [6,15,25]. We illustrate the approach with a random treatment procedure for infected nodes with degree k, such that they recover with an increased rate α+γwk, where γ is the overall treatment rate, and wk is a targeting fraction of the infected population with degree k: kwk=1.

In the weak limit we expect optimal treatment to favor large k (similar to targeted immunization), since yk* and pk* are proportional to k[26]. However, in the strong limit, treating low-degree nodes close to k will tend to decrease S, since their numbers must be lowered in the first step, before momenta differ significantly from zero. In intermediate cases, we expect a mixed strategy to minimize S. Treatment results are shown in Fig. 3 for a simple bimodal network with two degree classes for clarity, as a function of the targeting fraction for high-degree nodes, w50.

Action versus the fraction of infected high-degree nodes treated in a bimodal network [see Figs. 1(a) and 1(b)] and increasing treatment rate, γ (red→black): β˜=0.275. The inset shows the extinction times for a CMN with N=200.
FIG. 3.
Action versus the fraction of infected high-degree nodes treated in a bimodal network [see Figs. 1(a) and 1(b)] and increasing treatment rate, γ (redblack): β˜=0.275. The inset shows the extinction times for a CMN with N=200.

Interestingly, we find that choosing the optimal wk for the bimodal network can result in a nearly 50% decrease in the network action, implying an enormous reduction in extinction times, i.e., TT1/2 (inset of Fig. 3). Furthermore, we note that in Fig. 3 the size of the endemic state is minimized when w50=0 (circles), the equilibrium treatment rate γkwkgkyk*, is maximized when w50=0 (circles), and R0 is minimized when w50=1 (triangles), but none corresponds to the minimum extinction time control (squares) [27]. The example demonstrates that designing optimal controls intended to drive epidemics to extinction in finite networks cannot be found from the intuitive results of the deterministic limit alone, pk=0, but by targeting the network’s components in such a way as to minimize the network’s action.

In conclusion, we have considered how fluctuations in the SIS model produce extinctions from internal noise in finite heterogeneous networks, and have found that the process is captured by a most probable path. We were able to construct paths by combining the theory of rare events and random networks with a general degree distribution, and predict important consequences, such as the exponential decrease in extinction times with topological variation, as well as the multistep scaling of extinction through nodes with very different degree. Furthermore, we demonstrated how the theory can be used to manipulate fluctuations for optimal network control, producing exponential decreases in extinction times with a simple treatment strategy that minimized the action by leveraging its dependence on topology. Our theoretically predicted results were confirmed by simulations over large parameter ranges and different network topologies.

Lastly, we suggest the theoretical approach can be tailored to more arbitrary weighted networks and general epidemic processes that would allow one to predict the paths to extinction through real networks [1,8,13,17]. The specific formalism presented here could be augmented, in addition, to include degree correlations that may amplify or reverse the patterns described in interesting ways depending on the network assortativity.


We are grateful to Luis Mier-y-Teran Romero, D. J. Schneider, B. S. Lindley, C. R. Myers, and L. B. Shaw for useful discussions. J. H. received support as a National Research Council postdoctoral fellow. I. B. S was supported by U.S. Naval Research Laboratory funding (Grant No. N0001414WX00023) and the Office of Naval Research (Grant No. N0001414WX20610).



R. M. Anderson and R. M. May, Infectious Diseases of Humans (Oxford University Press, New York, 1991).


J. R. Banavar and A. Maritan, Nature (London)460, 334 (2009).NATUAS0028-0836


M. I. Dykman, I. B. Schwartz, and A. S. Landsman, Phys. Rev. Lett.101, 078101 (2008).PRLTAO0031-9007


M. Khasin and M. I. Dykman, Phys. Rev. Lett.103, 068101 (2009).PRLTAO0031-9007


A. Kamenev, B. Meerson, and B. Shklovskii, Phys. Rev. Lett.101, 268103 (2008).PRLTAO0031-9007


L. Billings, L. Mier-y Teran Romero, B. S. Lindley, and I. B. Schwartz, PLoS One8, e70211 (2013).POLNCL1932-6203


B. S. Lindley, L. B. Shaw, and I. B. Schwartz, Europhys. Lett.108, 58008 (2014).EULEEJ0295-5075


R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, Rev. Mod. Phys.87, 925 (2015).RMPHAT0034-6861


R. Pastor-Satorras and A. Vespignani, Phys. Rev. Lett.86, 3200 (2001).PRLTAO0031-9007


S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Rev. Mod. Phys.80, 1275 (2008).RMPHAT0034-6861


M. E. J. Newman, Phys. Rev. E66, 016128 (2002).PRESCM1063-651X


M. A. Muñoz, R. Juhász, C. Castellano, and G. Ódor, Phys. Rev. Lett.105, 128701 (2010).PRLTAO0031-9007


A. V. Goltsev, S. N. Dorogovtsev, J. G. Oliveira, and J. F. F. Mendes, Phys. Rev. Lett.109, 128702 (2012).PRLTAO0031-9007


E. Valdano, L. Ferreri, C. Poletto, and V. Colizza, Phys. Rev. X5, 021005 (2015).PRXHAE2160-3308


D. K. Wells, W. L. Kath, and A. E. Motter, Phys. Rev. X5, 031036 (2015).PRXHAE2160-3308


M. E. J. Newman, Networks: An Introduction (Oxford University Press, New York, 2010).


For example, analogous results can be derived by directly approximating A, in which case degrees and their statistical moments are replaced by eigenvector centralities and eigenvalues of A[8,13].


M. I. Dykman, E. Mori, J. Ross, and P. M. Hunt, J. Chem. Phys.100, 5735 (1994).JCPSA60021-9606


V. Elgart and A. Kamenev, Phys. Rev. E70, 041106 (2004).PRESCM1539-3755


M. I. Friedlin and A. D. Wentzell, Random Perturbations of Dynamical Systems, 2nd ed. (Springer-Verlag, New York, 1998).


B. S. Lindley and I. B. Schwartz, Physica (Amsterdam)255D, 22 (2013).PDNPDT0167-2789


C. R. Doering, K. V. Sargsyan, and L. M. Sander, Multiscale Model. Simul.3, 283 (2005).1540-3459, doi: 10.1137/030602800


P. Holme, PLoS One8, e84429 (2013).POLNCL1932-6203


Since the path is hyperbolic and power laws have many degree classes, it is useful to apply a binning procedure that reduces the dimension for the iterative action minimizing method. Above, we have uniformly binned the distribution kgk/k, and replaced k with the average in each bin. This was chosen because it gives accurate approximations for k2/k. Also, kmax was chosen small enough such that at least O(1) nodes could be expected in its bin given N.


K. Drakopoulos, A. Ozdaglar, and J. N. Tsitsiklis, IEEE. Trans. Netw. Sci. Eng.1, 67 (2014)., doi: 10.1109/TNSE.2015.2393291


R. Pastor-Satorras and A. Vespignani, Phys. Rev. E65, 036104 (2002).PRESCM1063-651X


With treatment, R0-1 is the distance from the critical point and satisfies k(βk2gk/k)/(R0+γwk/α)=1. Extinction and Control in Heterogeneous Networks&author=Jason Hindes,Ira B. Schwartz,&keyword=&subject=Letters,Polymer, Soft Matter, Biological, Climate, and Interdisciplinary Physics,