Rare events in networks with internal and external noise

Many complex systems of interest can display rare and extreme events in their dynamics. Though such events occur infrequently, their existence may have drastic consequences for a system's long-term function and behavior. Examples include recurrent dynamics in excitable media with noise [1], rogue waves in hydrodynamics and optics [2–4], fixation in evolutionary games [5], viral clearance [6], and switching in quantum mechanical oscillators [7]. For systems with explicit noise dependence, large-deviation theory (LDT) has been developed, which can be used to predict the most likely pathway through which a rare event occurs [8–10]. The approach approximates well a broad class of rare dynamical processes described in the exponential tail of a probability distribution [9,11]. Some recent work has begun to apply this formalism to high-dimensional problems in network prediction and optimal control [12–15]. On the other hand, there is interest in understanding the role of network topology on rare events in discrete stochastic processes, where noise is implicit and internal, and the direct application of LDT techniques is less straightforward [16–19]. In many situations, particularly in biology, a network may be subjected to both external environmental uncertainties and the inherent randomness of individual events [20–22]. Yet, no general methodology exists for analyzing rare processes in networks with internal and external noise sources – the problem addressed in this letter.

We first informally define a rare event as the appearance of some state in a network that is qualitatively different from its mean state, which occurs with exponentially small probability. Let us consider networks defined by a fixed set of (N) nodes and connections between them (called edges), such that they can be represented by adjacency matrices, A, whose elements are binary, i.e., $A_{ij}=1$ if nodes i and j are connected and zero otherwise. In the following, we study a simple stochastic process explicitly, which models infection spreading over a contact network: the Susceptible-Infected-Susceptible model (SIS) [23]. However, such assumptions are illustrative only and can be generalized to other dynamical processes (as explained below), weighted networks, adaptive networks, etc. See for instance [16,24,25].

In the SIS model nodes are divided into two types: infected and susceptible. Infected nodes transfer infection to their susceptible neighbors in a network at time, t, with probability-per-unit-time $\beta(t)$, and become susceptible with probability-per-unit-time $\alpha(t)$, where $\beta(t)$ and $\alpha(t)$ are known as the infection and recovery rates, respectively. In addition to the reactions on the network, $\beta(t)$ and $\alpha(t)$ may evolve according to independent stochastic processes. For illustration purposes, we take $\alpha(t)=\alpha=1$ and $\beta(t)=\beta_{0}[1-\xi(t)]$, where $\xi(t)$ is a zero-mean Gaussian white-noise process (GWNP) defined by the autocorrelation $\langle\xi(t)\xi(t')\rangle=2D\delta(t-t')$, and $\beta_{0}$ is the average infection rate. Fluctuations in $\beta(t)$ model temporal variations in environmental conditions, close-proximity contact rates, etc. [26,27].

As long as $\beta_{0}$ is above threshold, $\beta_{0}>\beta_{c}$, an endemic state emerges in large networks with a metastable number of infected nodes, after an initial transient period of seed infection and growth. The threshold $\beta_{c}$ can be calculated approximately with a variety of mean-field techniques [23]; below, we use the so called "pair quenched technique" (PQT), which is regarded as one of the most accurate for general network topologies [28].

Though infection is metastable above threshold, a network continually fluctuates in time due to random perturbations in the infection rate and changes in nodal states, which we call external and internal noise, respectively. Over time, noise drives the network to explore rare states, the most drastic of which is extinction (in the SIS model), where the number of infected nodes goes to zero [29,30]. Preliminary insight on how such rare states emerge can be gained by examining histograms of time-series data from stochastic realizations of extinction [25]. Example histograms are shown in fig. 1 for two networks with D = 0. As is typical for general networks, we find an exponential probability distribution of large deviations from the endemic state, which terminate at extinction. The distribution is peaked around the endemic state, and its shape depends on the infection rate and underlying network topology [16,30].

Zoom In Zoom Out Reset image size

A theory of rare events on networks should allow for predictions of the observed distribution of network states leading to extinction and the expected rate of extinction. Our analytical approach is to derive an approximate time-evolution equation for the probability of network states in the SIS model that is consistent with PQT, and analyze its maximum-probability solutions. We emphasize that much of the formalism developed below can be applied to other rare processes, such as network switching [24]. In fact, the primary difference among distinct classes of rare events resides in the boundary conditions needed to complete the theory and determine solutions of the LDT equations of motion (such as those derived below) [29,31].

As in [16,24], we introduce a network ensemble, for analysis purposes, consisting of C independent stochastic realizations of the SIS process on a fixed network. Ultimately, we are interested in the statistics of rare states in the limit of many realizations, or C → ∞. Following PQT [28], we describe the ensemble dynamics by the fraction of realizations in which node i is infected at time t, x_i(t), and the fraction of realizations in which node i is susceptible and node j is infected in an $i\text{-}j$ edge, $\phi_{ij}(t)$. These variables can be represented by an N-dimensional vector, $\vec{x}$, and an $N\times N$ dimensional matrix ϕ.

The probability distribution for $\vec{x}$ and ϕ in the ensemble at time t given $\beta(t)$, $P(\vec{x},\bm{\phi},t;\beta)$, is governed by a master equation (ME). The ME is a probability-flux equation that depends on the rates of all possible reactions in the ensemble and the increments to $\vec{x}$ and ϕ when a given reaction occurs. To see why the PQT is necessary, consider what happens when an infected node j infects a node i. The fractions x_i and $\phi_{ij}$ are incremented in a simple way: $x_{i}\rightarrow x_{i}+1/C$ and $\phi_{ij}\rightarrow\phi_{ij}-1/C$. However for additional neighbors of node i, e.g., node $l\neq j$, the increments to $\phi_{il}$ and $\phi_{li}$ depend on the state of l, which might be susceptible or infected. The probability that node l is infected is given by the fraction of realizations in the ensemble in which the triplet $l\text{-}i\text{-}j$ has an infected node l, a susceptible node i, and an infected node j. To keep our description limited to one and two-node variables, $\vec{x}$ and ϕ, we assume that the probabilities for states of nodes l and j in an $l\text{-}i\text{-}j$ triplet are independent.

The independence assumption enters when approximating the reaction increments as follows. When node j infects node i

$$\begin{align} \bm{\phi}\overset{i\text{-}j}{\underset{\text{infect}}\longrightarrow}\bm{\phi}+\bm{\mathcal{I}}^{ij}/C. \end{align} \tag{ 1 }$$

We approximate the increment to $\phi_{il}$, $\mathcal{I}^{ij}_{il}$, by minus the probability that node l is infected given node i is susceptible (independent of j). Similarly, we approximate the increment $\mathcal{I}^{ij}_{li}$ by the probability that node l is susceptible given node i is susceptible (independent of j):

$$\begin{eqnarray} \mathcal{I}^{ij}_{il}&= -\phi_{il}A_{il}A_{ij}(1-\delta_{jl})/(1-x_{i}), \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \mathcal{I}^{ij}_{li}&= (1-x_{i}-\phi_{il})A_{il}A_{ij}(1-\delta_{jl})/(1-x_{i}), \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} \mathcal{I}^{ij}_{ij}&= -A_{ij}, \quad \mathcal{I}^{ij}_{ji}=0. \end{eqnarray} \tag{ 4 }$$

On the other hand, if node i recovers, $x_{i}\rightarrow x_{i}-1/C$ and

$$\begin{align} \bm{\phi}\overset{i}{\underset{\text{recover}}\longrightarrow}\bm{\phi}+\bm{\mathcal{R}}^{i}/C. \end{align} \tag{ 5 }$$

We approximate the increment to $\phi_{ij}$, $\mathcal{R}^{i}_{ij}$, by the probability that node j is infected given node i is infected, and approximate the increment $\mathcal{R}^{i}_{ji}$ by minus the probability that node j is susceptible given node i is infected:

$$\begin{eqnarray} \mathcal{R}^{i}_{ij}&= (x_{i}-\phi_{ji})A_{ij}/x_{i}, \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} \mathcal{R}^{i}_{ji}&= -\phi_{ji}A_{ij}/x_{i}. \end{eqnarray} \tag{ 7 }$$

Altogether the increments give the ME for the network ensemble

$$\begin{eqnarray} & \frac{1}{C}\frac{\partial P}{\partial t}(\vec{x},\bm{\phi},t;\beta)\!=\!\alpha\sum_{i}\left[x_{i}\!+\!\frac{1}{C}\right]P\left(\vec{x}\!+\!\frac{\vec{1}_{i}}{C},\bm{\phi}\!-\!\frac{\bm{\mathcal{R}}^{i}}{C},t;\beta\right)\nonumber \\ & +\beta\sum_{ij}A_{ij}\left[\phi_{ij}+\frac{1}{C}\right]P\left(\vec{x}-\frac{\vec{1}_{i}}{C},\bm{\phi}-\frac{\bm{\mathcal{I}}^{ij}}{C},t;\beta\right)\nonumber \\ & -\alpha\sum_{i}x_{i}P\left(\vec{x},\bm{\phi},t;\beta\right)-\beta\sum_{ij}A_{ij}\phi_{ij}P\left(\vec{x},\bm{\phi},t;\beta\right), \end{eqnarray} \tag{ 8 }$$

where $\vec{x}+\vec{1}_{i}/C=(x_{1},x_{2},\ldots,x_{i}+1/C,\ldots,x_{N})$.

We would like to extract a solution of eq. (8) that describes rare events. The functional form for the probability of network states in individual stochastic realizations is suggested by fig. 1. As remarked above such distributions, $\rho(\vec{x},\bm{\phi},t;\beta)$, are exponential functions of $\vec{x}$ and ϕ,

$$\begin{equation} \rho(\vec{x},\bm{\phi},t;\beta)=b\exp\{-S(\vec{x},\bm{\phi},t;\beta)\}, \end{equation} \tag{ 9 }$$

where b is a normalization constant. As a consequence, the ensemble distribution takes a Wentzel-Kramers-Brillouin form (WKB), $P(\vec{x},\bm{\phi},t;\beta)=a\exp\{-CS(\vec{x},\bm{\phi},t;\beta)\}$: a product of independent and identical distributions. The WKB form is an essential property of the probability for rare processes describable by LDT [9]. Substituting the WKB ansatz into eq. (8), Taylor expanding eq. (8) in powers of the small parameter $1/C$, and neglecting terms of $\mathcal{O}(1/C)$ or smaller [16,24,32], gives a Hamilton-Jacobi equation for the probability exponent $S(\vec{x},\bm{\phi},t;\beta)$:

$$\begin{equation} \frac{\partial S}{\partial t} + H\Big(\vec{x},\bm{\phi},\frac{\partial S}{\partial \vec{x}},\frac{\partial S}{\partial \bm{\phi}};\beta\Big)=0. \end{equation} \tag{ 10 }$$

As in analytical mechanics, S and H are known as the action and Hamiltonian, respectively. In particular, the Hamiltonian is a function of derivatives of S, defined as conjugate momenta

$$\begin{equation} p_{i}=\frac{\partial S}{\partial x_{i}}, \;\;\text{and}\;\; m_{ij}=\frac{\partial S}{\partial\phi_{ij}}, \end{equation} \tag{ 11 }$$

with

$$\begin{eqnarray} & H(\vec{x},\bm{\phi},\vec{p},\bm{m};\beta)\!=\!\alpha\sum_{i}x_{i}\bigg(\exp\Bigg\{\!-\!p_{i}\!+\!\sum_{j}\frac{A_{ij}}{x_{i}}\bigg[\!-\!\phi_{ji}m_{ji}\nonumber \\ & \,+(x_{i}-\phi_{ji})m_{ij}\bigg]\Bigg\}-1\bigg)+\beta\sum_{ij}A_{ij}\phi_{ij}\bigg(\exp\Bigg\{p_{i}-m_{ij}\nonumber\\ & \,+\sum_{l}A_{il}\frac{(1-\delta_{jl})}{(1-x_{i})}\big[-\phi_{il}m_{il}+(1-x_{i}-\phi_{il})m_{li}\big]\Bigg\}-1\bigg). \end{eqnarray} \tag{ 12 }$$

It is important to note that the elements of ϕ are not independent. This can be seen by considering the fraction of realizations in the ensemble in which both nodes sharing an edge are infected. For example, the fraction of such $i\text{-}j$ edges is $x_{j}-\phi_{ij}=x_{i}-\phi_{ji}$ [28], which implies $\phi_{ji}=x_{i}-x_{j}+\phi_{ij}$. If such a constraint is enforced from the onset, it implies that $m_{ji}=0$ in eq. (12), since the probability does not depend on $\phi_{ji}$ explicitly. We emphasize that when constructing a Hamiltonian from eq. (12) using such constraints in ϕ, one must select a permutation for each $i\text{-}j$ edge in the network, e.g., $\phi_{ij}$, replace the other permutation $\phi_{ji}$ by $x_{i}-x_{j}+\phi_{ij}$, and set $m_{ji}=0$ in eq. (12). This procedure is useful for eliminating redundancy. We demonstrate with examples shortly.

The solution of eq. (10) for the probability exponent (action) is well known from analytical mechanics,

$$\begin{equation} S(\vec{x},\bm{\phi},t;\beta)=\int\big[\vec{p}\cdot\dot{\vec{x}}\;+\langle\bm{m},\dot{\bm{\phi}}\rangle-H(\vec{x},\bm{\phi},\vec{p},\bm{m};\beta)\big]\text{d}t, \end{equation} \tag{ 13 }$$

with scalar products $\vec{p}\cdot\dot{\vec{x}}=\sum_{i}p_{i}\dot{x}_{i}$ and $\langle\bm{m},\dot{\bm{\phi}}\rangle=\sum_{ij}m_{ij}\dot{\phi}_{ij}$, and where $\vec{x}$, $\vec{p}$, ϕ, and m satisfy Hamilton's equations of motion [29,32]:

$$\begin{equation} \dot{x_{i}}=\frac{\partial H}{\partial p_{i}}, \;\; \dot{\phi}_{ij}=\frac{\partial H}{\partial m_{ij}}, \;\; \dot{p_{i}}=-\frac{\partial H}{\partial x_{i}}, \;\; \dot{m}_{ij}=-\frac{\partial H}{\partial \phi_{ij}}. \end{equation} \tag{ 14 }$$

In general, rare events are characterized by ($\vec{x}$, $\vec{p}$, ϕ, m)-trajectories of the LDT equations of motion (e.g., eq. (14)) with particular boundary conditions, and their probabilities given by eq. (9) and eq. (13) to logarithmic accuracy. Since such solutions minimize the action, they are most probable, or optimal solutions [32].

In the case of extinction from a metastable endemic state with constant infection rate, $\beta=\beta_{0}$, the solution has a special form, since there is approximately no time-dependence to $\rho(\vec{x},\bm{\phi};\beta)$ [30]. Hence $\partial S/\partial t=H(\vec{x},\bm{\phi},\vec{p},\bm{m};\beta)=0$ [20,32]. Metastability, therefore, implies that network extinction is a zero-energy solution of eq. (14) that starts at an endemic state (a fixed point of eq. (14) with $\vec{x}\neq\vec{0}$, $\bm{\phi}\neq\bm{0}$) and ends at the extinct state (a fixed point of eq. (14) with $\vec{x}=\vec{0}$, $\bm{\phi}=\bm{0}$). These boundary conditions are sufficient to determine solutions numerically [33]² .

In order to highlight features of the combined WKB and PQT techniques, we look at two example networks where $\beta=\beta_{0}$, before moving on to the case where external noise is present. The first example is a random network where all nodes have k neighbors (k-regular) [28]. Such networks are useful models of homogeneous topologies with local interactions, and where the fully connected graph can be recovered as $k\rightarrow N-1$. As written, the dynamical system eq. (14) is high dimensional. However, since all nodes have the same degree, we can significantly reduce the dimension by approximating local variables by their averages: $x=\sum_{i}x_{i}/N$, $\phi=\sum_{ij}\phi_{ij}/[Nk]$, $p=\sum_{i}p_{i}/N$, and $m=\sum_{ij}m_{ij}/[Nk]$. Assuming that $x_{i}=x$, $\phi_{ij}=\phi$, $p_{i}=p$, $m_{ij}=m$ for all i and j in the network, the Hamiltonian becomes $H(\vec{x},\bm{\phi},\vec{p},\bm{m};\beta)=H_{k}(x,\phi,p,m;\beta)$ with

$$\begin{eqnarray} & H_{k}(x,\phi,p,m;\beta)=N\alpha x\left[\exp\left\{-p+km(x-2\phi)/x\right\}-1\right]\nonumber\\ & +Nk\beta\phi[\exp\{p-m+\!(k-1)m(1-x-2\phi)/(1-x)\}\!-\!1]. \end{eqnarray} \tag{ 15 }$$

From eq. (14) it follows that $\dot{x}=N^{-1}\partial H/\partial p$, $\dot{p}= -N^{-1}\partial H/\partial x$, $\dot{\phi}=[Nk]^{-1}\partial H/\partial m$, and $\dot{m}=-[Nk]^{-1}\partial H/\partial \phi$. Hence, the equations of motion for $\dot{x}$ and $\dot{\phi}$ are:

$$\begin{eqnarray} \dot{x}&= -\alpha x \exp\big\{-p+km(x-2\phi)/x\big\}+\beta k\phi\exp\big\{p-m\nonumber \\ &\quad +\,(k-1)m(1-x-2\phi)/(1-x)\big\}, \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} \dot{\phi}&= \alpha\big[x-2\phi\big]\exp\big\{-p+km(x-2\phi)/x\big\}\nonumber \\ &\quad +\,\beta\phi\left[\frac{(k-1)(1-x-2\phi)}{(1-x)}-1\right]\exp\big\{p-m\nonumber\\ &\quad +\,(k-1)m(1-x-2\phi)/(1-x)\big\}. \end{eqnarray} \tag{ 17 }$$

We point out that the known pair-approximation equations for k-regular networks are found by setting the momenta equal to zero in eqs. (16), (17), $p=m=0$ [28]. In fact, this is a general feature of our approach: network mean-field theories are zero-momentum invariant manifolds of Hamilton's equations of motion [24]. In general, rare events are describable as trajectories that lie off of this manifold, i.e., $\vec{p}\neq\vec{0}$ and $\bm{m}\neq\bm{0}$ [29]³.

The second example is a star network ($\mathcal{S}(L)$), which is composed of a single central node sharing an edge with $L=N-1$ "leaf" nodes. Leaf nodes do not share edges. Star networks are relevant because they exhibit localization of the endemic state around the central node near threshold. Localization is ubiquitous and is observed, for instance, in certain classes of random networks with unbounded degree distributions as N → ∞ [30,34–36].

Let us denote the infection and momentum of the central node with index 0, x₀ and p₀, respectively. Similar to the k-regular networks, we define average coordinates for the leaves: $x_{L}=\sum_{j\in\mathcal{L}}x_{j}/L$, $\phi_{L0}=\sum_{j\in\mathcal{L}}\phi_{j0}/L$, $\phi_{0L}={\sum_{j\in\mathcal{L}}\phi_{0j}/L}$, $p_{L}=\sum_{j\in\mathcal{L}}p_{j}/L$, $m_{L0}=\sum_{j\in\mathcal{L}}m_{j0}/L$, and $m_{0L}=\sum_{j\in\mathcal{L}}m_{0j}/L$, where $\mathcal{L}$ denotes the set of leaves. Since there are two types of nodes and only one type of edge in a $\mathcal{S}(L)$–between the central node and each leaf–we can apply the constraint $\phi_{0L}=x_{L}-x_{0}+\phi_{L0}$ with $m_{0L}=0$ (see footnote ³). Approximating all leaf variables by their respective averages gives the $\mathcal{S}(L)$ Hamiltonian

$$\begin{eqnarray} & H_{\text{star}}(\vec{x},\bm{\phi},\vec{p},\bm{m};\beta)=\alpha x_{0}\big[\exp\{-p_{0}-L\phi_{L0}m_{L0}/x_{0}\}-1\big]\nonumber\\ & +\,\alpha Lx_{L}\big[\exp\{-p_{L}+(x_{0}-\phi_{L0})m_{L0}/x_{L}\}-1\big]\nonumber \\ & +\,\beta L(x_{L}\!-\!x_{0}+\phi_{L0})\big[\exp\{p_{0}\!+\!(L\!-\!1)(1\!-\!x_{L}\!-\!\phi_{L0})m_{L0}\nonumber\\ & \,/(1-x_{0})\}-1\big]+\beta L\phi_{L0}\big[\exp\{p_{L}-m_{L0}\}-1\big]. \end{eqnarray} \tag{ 18 }$$

Predictions of $\rho(\vec{x},\bm{\phi},t;\beta)$ from solving eq. (14) for both networks agree well with measured histograms, fig. 1 (blue). In addition, fig. 2 shows path projections of many Monte Carlo simulations that precede extinction for both networks (see footnote ²). Panel (a) is a heat map for a 6-regular network, giving the extinction-path density as a function of the fraction of infected nodes and edges. The predicted path is in blue, which lies on the maximum of the density, as expected. The dashed green line represents a path that assumes statistical independence of neighbors, $\phi=x(1-x)$, which is consistent with standard heterogeneous and quenched mean-field approximations [16,24]; such a path does not lie on the maximum density. Therefore, we see that the PQT provides a significant quantitative improvement over standard mean-field approximations in predicting extinction [28].

Zoom In Zoom Out Reset image size

Panel (b) shows the fraction of 10⁴ stochastic simulations in which the central node of a $\mathcal{S}(L)$ is infected while a fraction x_L leaf nodes are infected, just before extinction occurs (red circles). The predicted path from eq. (14) and eq. (18) is shown in blue, and agrees qualitatively well with simulation results. The red-dashed curve is the predicted path into the endemic state $(p_{0}=p_{L}=m_{L0}=0)$ from an initial seed infection, $x_{L}=\epsilon\ll1$. The dissimilarity between the two paths (only overlapping at the endpoints) demonstrates the highly non-reversible nature of rare events in heterogeneous networks driven by internal noise [14,16].

In addition to the paths and probabilities of rare events, we are also interested in their frequency (rate). Generally, rare events are Poisson processes with rates proportional to their probabilities [9,20,37]. Hence the inverse rate of extinction, or average time $\langle T\rangle$, is expected to be inversely proportional to the probability of extinction, or to logarithmic accuracy

$$\begin{equation} \ln\langle T\rangle=S(\vec{x}=\vec{0},\bm{\phi}=\bm{0};\beta)+B, \end{equation} \tag{ 19 }$$

where B is assumed to be an $\mathcal{O}(1)$ pre-factor [31,37]. Figure 3 shows the average extinction time as a function of $\beta/\beta_{cr}$ [28] for a 6-regular network and a $\mathcal{S}(L)$ with D = 0. Predictions from eq. (19) are shown in red with good agreement.

Zoom In Zoom Out Reset image size

So, far we have considered only the effects of internal noise on the occurrence of rare extinction. If $\beta(t)$ has a non-zero noise amplitude, then we must integrate over realizations of $\beta(t)$ to get the marginal distribution of rare network states, $\rho(\vec{x},\bm{\phi},t)$. In the case where $\beta(t)=\beta_{0}[1-\xi(t)]$, and $\xi(t)$ is a GWNP, the statistical weight for a given realization of $\xi(t)$ is $\mathcal{P}[\xi(t)]\sim\exp\{-\int[\xi^{2}/4D]\text{d}t\}$ [20]. Therefore, the marginal distribution takes the functional form:

$$\begin{eqnarray} \rho(\vec{x},\bm{\phi},t)&= \int\mathcal{D}\xi\;\mathcal{P}[\xi(t)]\rho(\vec{x},\bm{\phi},t;\beta(\xi))\nonumber\\ &= \int\mathcal{D}\xi\exp\left\{\!-\!\int\left[\vec{p}\cdot\dot{\vec{x}}\;\!+\!\langle\bm{m},\dot{\bm{\phi}}\rangle\right.\right.\nonumber\\ & \quad-\left.\left.H(\vec{x},\bm{\phi},\vec{p},\bm{m};\beta)\!+\!\frac{\xi^{2}}{4D}\right]\text{d}t\right\}. \end{eqnarray} \tag{ 20 }$$

As we have seen, rare events occur with exponentially small probabilities along optimal paths described by analytical mechanics. Leveraging this insight once again, we approximate eq. (20) by its maximum-probability contribution. We note that the exponent in eq. (20) takes the form of an effective action, $S_{\text{eff}}(\vec{x},\bm{\phi},\xi,t)$, whose effective Hamiltonian is

$$\begin{equation} H_{\text{eff}}(\vec{x},\bm{\phi},\vec{p},\bm{m},\xi)=H(\vec{x},\bm{\phi},\vec{p},\bm{m};\beta(\xi))-\frac{\xi^2}{4D}. \end{equation} \tag{ 21 }$$

Therefore, the minimum action solution satisfies another Hamilton-Jacobi equation, $\partial S_{\text{eff}}/\partial t +H_{\text{eff}}=0$, with ${H_{\text{eff}}=0}$.

When external noise is a GWNP and $H_{\text{eff}}(\vec{x},\bm{\phi},\vec{p},\bm{m},\xi)$ is quadratic in ξ, we can solve explicitly for the realization of ξ with maximum probability as a function of the optimal path coordinates, $\xi(t)=\bar{\xi}(\vec{x},\bm{\phi},\vec{p},\bm{m})$. Substituting, $\beta=\bar{\beta}=\beta_{0}[1-\bar{\xi}]$ into eq. (12) and eq. (21), with the latter set equal to zero, we find⁴:

$$\begin{eqnarray} & \,\bar{\xi}=-2D\beta_{0}\sum_{ij}A_{ij}\phi_{ij}\bigg(\exp\Big\{p_{i}-m_{ij} \\ & \,+\sum_{l}A_{il}\frac{(1-\delta_{jl})}{(1-x_{i})}\big[-\phi_{il}m_{il}+(1-x_{i}-\phi_{il})m_{li}\big]\Big\}-1\bigg).\nonumber \end{eqnarray} \tag{ 22 }$$

As with the D = 0 case, optimal paths with external noise satisfy eq. (14). The main difference is that the infection rate implicit in eq. (14) is set equal to the optimal realization, $\beta(t)=\bar{\beta}$ (see footnote ⁵). The fixed-point boundary conditions unchanged. Just as before, the action is computed along the path,

$$\begin{equation} S_{\text{eff}}\big(\vec{x},\bm{\phi};\bar{\beta}\;\big)=\int\big[\vec{p}\cdot\dot{\vec{x}}\;+\langle\bm{m},\dot{\bm{\phi}}\rangle\big]\text{d}t, \end{equation} \tag{ 23 }$$

giving average extinction times:

$$\begin{equation} \ln\langle T\rangle=S_{\text{eff}}(\vec{x}=\vec{0},\bm{\phi}=\bm{0};\bar{\beta})+B. \end{equation} \tag{ 24 }$$

Predictions of $\ln\langle T \rangle$ from eq. (24) are compared to simulations in fig. 4, with changing D and fixed $\beta_{0}$, and show good agreement for the two example networks (see footnote ²). Qualitatively, we find that $\ln\langle T\rangle$ decreases as D increases —demonstrating that large fluctuations producing extinctions in networks are only enhanced by external noise. In addition, the shape of the optimal external noise in simulations can be compared directly to predictions, eq. (22). For example, fig. 5 shows the generic form of ξ, in which the external noise increases from zero at the endemic state, reaches a maximum at some intermediate infection level, and then decreases to zero at the extinct state.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given the outlined mechanics for rare dynamics in networks with internal and external noise, we may wonder how each differs in its effect on extinction. For instance, it is well known that the probability of large deviations from internal noise generally decreases exponentially with the network size, N, such that $S\sim N$ as N → ∞ [9,16,32]. An interesting question concerns how $S_{\text{eff}}(\vec{x}=\vec{0},\bm{\phi}=\bm{0};\bar{\beta})$ (or $S_{\text{eff}}$ for short) scales with N in the presence of external noise and different network topology. Figure 6 shows $S_{\text{eff}}$ as a function of N for two example networks and several values of D. Panel (a) shows results for an 8-regular network. For $D\neq0$, we can see that $S_{\text{eff}}$ scales linearly with N (black dashed line) for small N, implying that for sufficiently small networks, internal noise is dominant and external noise has a only a perturbative effect. Therefore for k-regular networks there is always some range of N over which $S_{\text{eff}}$ scales linearly with N for any D. As N increases, however, their is a crossover to $S_{\text{eff}}\sim1/D$, as the intensity of internal fluctuations becomes negligible and extinction depends effectively on fluctuations in β alone [22].

Zoom In Zoom Out Reset image size

On the other hand, panel (b) shows the result for a $\mathcal{S}(L)$. In contrast, there is a minimum network size for fixed $\beta_{0}$, $N_{c}(\beta_{0})$, such that smaller stars cannot stabilize endemic states. As a consequence, $S_{\text{eff}}$ grows nonlinearly from zero for $N \gtrsim N_{c}$ (black dashed line), even when D = 0. The linear scaling (red dashed line) is found when D = 0, but only asymptotically as N → ∞, as the dashed lines approach each other. The internal noise scaling pattern for a $\mathcal{S}(L)$ is a consequence of localization near threshold, and we expect it to appear in networks that show localization more generally (see footnote ³). For $D\neq0$, there is a crossover from the internal noise scaling pattern to $S_{\text{eff}}\sim1/D$ as N → ∞, just as in the k-regular case. However, since the linear scaling with N is reached only asymptotically in the absence of external noise, it is possible that such scaling is never found if the external noise is too large, e.g., $D=0.05$ in fig. 6(b).

We conjecture that the noise crossover for localized endemic states should be observed in large power-law networks with degree exponents greater than three. This is because the endemic state near threshold⁶ in such networks is known to behave like a star sub-graph of the full network with $L\rightarrow k_{max}$, where k_max is the maximum degree of the hub node [28].

In conclusion, we have developed an analytical framework that describes the most-probable pathway of rare events, such as extinction, in networks having both internal and external noise. Our formalism combined pair-quenched mean-field approximations for network dynamics with WKB techniques and large deviation theory. Using these techniques, we were able to predict the distribution of rare network states and the frequency of extinction, depending on both a network's topology and the intensity of external noise. We showed that the presence of any external noise leads to a crossover region where the familiar exponential dependence of large deviations on a network's size (N) is lost —in the infinite-size limit, depending only on the external noise intensity, D.

We observed that the crossover occurs as a network's probability exponent approaches a scale of $1/D$. As a consequence, there may be no range of N over which the exponential scaling with N is found for topologies that are localized near threshold, since such networks show the exponential scaling with N only asymptotically as N → ∞ (without external noise). With external noise, the crossover may occur for N that is too small to show exponential scaling.

JH is a National Research Council postdoctoral fellow. IBS was supported by the U.S. Naval Research Laboratory funding (N0001414WX00023) and office of Naval Research (N0001416WX00657) and (N0001416WX01643). We thank L. Mier-Y-Teran-Romero for useful discussions.