Epidemic Extinction and Control in Heterogeneous Networks.

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.

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 , where the degree is the number of links of a node. Simple graphs with nodes can be generated from in several ways, for example, with a configuration model network (CMN) [16]. Such networks are usefully represented by an adjacency matrix , where is 1 if nodes and 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 is either infected, denoted , or susceptible, , and changes its state with probability per unit time , and with probability per unit time , 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 with its expectation value in an ensemble of networks , which for uncorrelated networks takes the form in the limit of large . 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 , , which has corresponding reactions and rates: with rate of infection , and with recovery rate , where , and .

Since the SIS model is stochastic, associated with all network states is a probability distribution , which for networks having a general degree distribution satisfies an approximate master equation where . We are interested in the behavior of Eq. (1) for large, but finite, networks. As customary, we assume is large and take the leading order in a expansion [3,18]. This is similar to the Wentzel-Kramers-Brillouin (WKB) ansatz of quantum mechanics, where plays the role of Planck’s constant in Schrödinger’s equation. In accordance with WKB, by writing , where , and taking the leading order in , or and , we find a Hamilton-Jacobi equation, , where and are called the action and Hamiltonian, respectively. As in classical mechanics, the Hamiltonian is a function of the coordinate and its conjugate momentum : In this context, momenta behave as fluctuations on —describing both size and direction.

It is convenient to analyze Eq. (2) by solving the canonical equations of motion, , , in terms of the fraction of each degree class infected , the ratio , and the rescaled time :

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 is quasistationary, or in the WKB ansatz, : 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, , , to extinction with nonzero momentum, , [3,5,19]. The functions and depend on , with as , and 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].

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 , , and , and with two degree classes, ; nodes occupying 10% of the network. (b) A corresponding CMN. (c) A CMN with , , and where high-degree nodes have , and low-degree have . (d) A CMN with , , and , where high-degree nodes have , and low-degree have .

In general, the average time to reach extinction will depend on and network properties in complicated ways [22,23]. However, for sufficiently large , the transition is an exponential process with a rate proportional to the probability, and therefore . 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 (which we call “weak”), the paths to extinction are approximately linear from to . The explicit form can be seen by expanding the equilibria in powers of , which gives to first order , 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, (), 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 and , respectively, and in Fig. 2(a) (red).

FIG. 2.

(a) Projections of the optimal paths for the truncated power law [see Fig. 1(d)] shown for increasing (, ) in steps of 0.5, compared with the path into the endemic state for (green). Arrows indicate direction in time. (b) Projections into for the same distribution with and 17 bins [24], shown for bins with increasing (), 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 [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, (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 to linear order and studying the asymptotic scaling of for large . Inserting these assumptions into Eq. (3), we find that must satisfy an eigenvalue equation for the rate :

Since and the sum are independent, it must be that for large , we have [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 -dependent contribution to approaches a maximum along the path. In analogy with mechanics, this can be thought of as the network’s contribution to the “force” on , which near its maximum quickly “pushes” from near zero to its maximum . On the other hand, since and decreases along the path, by inspecting , we can find an upper bound for , , because the -dependent contribution to is positive. Since until differs from zero, as the force approaches its maximum, can be approximated by the upper bound, giving the scaling for large : .

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 , which gives an eigenvalue equation for the rate which in the limit of large implies . Examples are shown in Figs. 2(a) (black) and 2(b) (black); in the former, the strong limit scaling starts to be visible for . 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 , such that they recover with an increased rate , where is the overall treatment rate, and is a targeting fraction of the infected population with degree : .

In the weak limit we expect optimal treatment to favor large (similar to targeted immunization), since and are proportional to [26]. However, in the strong limit, treating low-degree nodes close to will tend to decrease , 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 . 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, .

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, (): . The inset shows the extinction times for a CMN with .

Interestingly, we find that choosing the optimal for the bimodal network can result in a nearly 50% decrease in the network action, implying an enormous reduction in extinction times, i.e., (inset of Fig. 3). Furthermore, we note that in Fig. 3 the size of the endemic state is minimized when (circles), the equilibrium treatment rate , is maximized when (circles), and is minimized when (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, , 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.