Stochasticity of infectious outbreaks and consequences for optimal interventions

The most commonly used model in epidemiology is the classical SIR model [10, 11]. It has been used to study the behavior of epidemics because it captures most of the general features of the dynamics and requires few parameters. The SIR model consists of a population of individuals, each of which can be susceptible (S), infected (I) or recovered (R). In the following, we use an extension of the SIR model, the so-called SEIR model, which contains an additional compartment of exposed (E) individuals. The exposed compartment contains infected individuals which are not yet able to transmit the infection, because the pathogen is in an incubation or pre-contagious phase. The SEIR model has previously been used to model epidemics on contact networks in the context of the 2005 SARS-CoV outbreak [27] and to estimate the turning period of SARS-CoV-2 outbreaks [28].

The deterministic SEIR dynamics is governed by a set of differential equations,

$$\begin{equation}\dot{S}=-\beta \frac{SI}{N}\end{equation} \tag{ 1 }$$

$$\begin{equation}\dot{E}=\beta \frac{SI}{N}-\sigma E\end{equation} \tag{ 2 }$$

$$\begin{equation}\dot{I}=\sigma E-\gamma I\end{equation} \tag{ 3 }$$

$$\begin{equation}\dot{R}=\gamma I.\end{equation} \tag{ 4 }$$

These equations describe the infection dynamics in a well-mixed population of constant size N. Susceptible individuals get exposed at a rate proportional to the fraction of infected individuals with a transmission rate parameter β. Exposed individuals become infectious at a constant rate σ; this delay represents intra-host growth of the viral population up to a level where transmission to a next host becomes likely. Infected individuals recover at a constant rate γ (figure 1(a)). In the model tuned to SARS-CoV-2, we use estimates of the average pre-contagious period (τ_σ ≡ σ⁻¹ = 4 d) and infectious period (τ_γ ≡ γ⁻¹ = 6 d), and we choose the parameter β from the interval 0.075–0.75 d⁻¹ [1, 29–31]. A further, frequently used parameter to characterize epidemic growth is the so-called reproductive number, which is defined a the average number of transmissions from a single individual during its infectious period. In the well-mixed SEIR model, the reproductive number is simply related to the basic rate parameters, R₀ = β/γ [17]; however, this relation is no longer valid on a contact network.

Zoom In Zoom Out Reset image size

We are interested in the first stage of the dynamics, where the approximation S ≈ N is valid and the number of infected individuals grows exponentially with rate λ, as given by

$$\begin{equation}\lambda ({R}_{0},\sigma ,\gamma )=-\frac{\sigma +\gamma }{2}+\frac{1}{2}\sqrt{{(\sigma -\gamma )}^{2}+4\sigma \gamma {R}_{0}}\end{equation} \tag{ 5 }$$

(see methods). Here we have used R₀ instead of β as the independent parameter, because the relation (5) remains valid on a contact network (see below). The reproductive number delineates the regimes of deterministic growth (λ > 0 for R > 1) and decline (λ < 0 for R₀ < 1).

In figure 1(b), we show the deterministic growth of an SEIR epidemic starting from a single infected individual. After a brief initial period, the infected and exposed populations grow at the same rate λ and at comparable relative size. Because incubation slows down growth, this rate is lower than for a SIR model at the same value of R₀ = β/γ, which is given by λ(R₀, γ) = (R₀ − 1)γ (dashed lines in figure 1(b)). In other words, in the deterministic dynamics, a sizeable exposed phase has a deleterious fitness effect for the pathogen. In the stochastic theory, as we discuss below, the exposed phase also generates a beneficial effect, because it reduces fluctuations leading to early extinctions.

Given initially low numbers of exposed and infected individuals, stochasticity plays an important role in the transmission dynamics of local outbreaks. In figure 2, we plot time courses of the number of infected individuals from simulations of the SEIR model in a well-mixed population. Some of these time courses reach only small transient population sizes and then go extinct (red lines), the remainder proliferate to deterministic, initially exponential growth (blue lines). Three characteristics of these stochastic dynamics are readily recognized. First, the fraction, duration, and population size of extinction-bound time courses strongly depends on the reproductive number. Below, we calculate the proliferation probability Q_p(R₀), given the initial condition of a single infected individual. Second, the fraction of extinction-bound trajectories decreases with increasing population size. Beyond a threshold n*(R₀), the so-called establishment size (black dashed lines in figure 2), proliferation becomes the more likely fate. The establishment size will also be calculated from a suitably extended expression of the extinction probability. Third, in outbreaks bound for proliferation, the average number of infected individuals (thick blue lines) grows initially faster than expected from the deterministic dynamics (blue dashed lines). This initial boost is a posterior effect of sampling only time courses leading to epidemic growth; it is more pronounced for values of R₀ close to 1 (figure 2(b)). In this regime, typical stochastic trajectories reach the epidemic regime (black dashed lines) several days earlier than expected from the deterministic model. These stochastic initial-time effects can drastically affect the efficiency of monitoring and control strategies.

Zoom In Zoom Out Reset image size

To study the stochastic dynamics of the SEIR model we assume that, at a given time, each individual j in the populations is in a state x_n ∈ {s, e, i, r}. The state of an individual at time t + dt, given the state at time t, is governed by the transition probabilities

$$\begin{equation}{\mathcal{P}}_{j\,}(s\to e)=\beta \frac{I}{N}\,\mathrm{d}t\end{equation} \tag{ 6 }$$

$$\begin{equation}{\mathcal{P}}_{j\,}(e\to i)=\sigma \,\mathrm{d}t\end{equation} \tag{ 7 }$$

$$\begin{equation}{\mathcal{P}}_{j\ }(i\to r)=\gamma \,\mathrm{d}t.\end{equation} \tag{ 8 }$$

In a well-mixed population, the transition rates are uniform across individuals and equal the deterministic rates in equations (1)–(4).

The proliferation probability of the well-mixed SEIR model can readily be computed from a branching process approximation with constant linear birth and death rates of infected individuals, which is appropriate in the exponential growth regime (S ≈ N). The well-known result for the SIR model, 1 − Q_p = 1/R₀ carries over to the SEIR model, because the probability of infections is independent of the incubation step [17, 18].

To compute the establishment size, we evaluate the conditional proliferation probability Q_p(R₀, n), given that the infected population has already reached a size n. This probability is a decreasing function of n, interpolating between the initial value Q_p(R₀, n) = Q_p(R₀) and near-certain epidemic growth, Q_p(R₀, n) ≃ 1, of sufficiently large infected clusters (figure 2). We obtain ${Q}_{\text{p}}({R}_{0},n)=1-{\left(1/{R}_{0}\right)}^{n}$, assuming that extinction from an initial state of n infected individuals is equivalent to n independent extinction processes, each of them starting from a single individual. The establishment size n*(R₀) is then obtained from the condition Q_p = 1/2, which gives n*(R₀) = log 2/log R₀ (black dashed lines in figures 2 and 4(b) and 6(b)).

Multiple stochastic factors govern the initial stage of an outbreak [13, 32–34]. In addition to the intrinsic stochasticity of the dynamics given by equations (6)–(8), there is heterogeneity in the initial condition of patient zero. Most studies have addressed differences in the reproductive number R₀ across the population, which can be caused by a complex mixture of host, pathogen and environmental factors [32]. Here, we consider stochastic effects generated by a network of contacts that controls the interactions between individuals and translates into an effective variation of infectiousness (figure A1(b)). This (undirected) network is characterized by a broad degree distribution $\hat{p}(k)\sim {k}^{-\alpha }$ [12], which gives the probability that an individual in the network has k contacts that can become transmission channels. While the degree distribution ignores the specific topology of the network, it captures the effect of largely connected nodes or hubs on the transmission statistics. Hubs play an important role in super-spreading events, which turns out to be the most relevant network feature for the degree-based testing protocol that we discuss below.

Given that infections take place between neighboring nodes of the contacts of the contact network, the SEIR transition probability from susceptible to exposed is given by

$$\begin{equation}{\mathcal{P}}_{j\,}(s\to i)=\beta \frac{{I}_{j}}{\langle k\rangle }\,\mathrm{d}t,\end{equation} \tag{ 9 }$$

where I_j is the number of infected contacts (neighboring nodes) of node j and $\left\langle \cdots \,\right\rangle$ denotes averages over the degree distribution $\hat{p\,}(k)$. This local process replaces the infection step of equation (7) in a well-mixed population. The normalization factor ⟨k⟩ in equation (9) ensures that in the limit of large ⟨k⟩ at constant β, the network dynamics converges to the well-mixed case.

In figure 3, we compare the SEIR infection dynamics on a contact network with the dynamics in a well-mixed population at the same parameters β, γ, σ. In the parameter regime shown, the epidemic on a network grows faster; most strikingly, the regime of positive growth extends well below the growth threshold β/γ = 1 in the well-mixed case (figure 3(a)). Consistently, the proliferation probability Q_p remains positive in the same regime; for larger values of β/γ, however, proliferation is less likely on a network than in a well-mixed population. This pattern indicates two opposing effects of the contact network on epidemic growth. First, the epidemic spreads preferentially on high-connectivity nodes, which enhances the effective transmission rate. Second, infections deplete susceptible individuals specifically in the neighborhood of an infected node, leading to saturation effects at larger reproductive numbers.

Zoom In Zoom Out Reset image size

Most importantly, the proliferation probability of outbreaks on a contact network is no longer uniform but depends strongly on the contacts of the first infected individual (patient zero). In figure 4(a), we plot Q_p as a function of the degree of patient zero, k₀, and the reproductive number. Similar to the population threshold defined before, we define a threshold degree ${k}_{0}^{\ast }({R}_{0})$, such that ${Q}_{\text{p}}({R}_{0},{k}_{0}^{\ast })=1/2$ (dashed line). As expected, highly connected seeding nodes in the network are more likely to produce epidemics than nodes with low connectivity. This is the so-called super-spreader effect: most of the epidemics are seeded by highly connected individuals. The effect is most pronounced for small values of R₀. Even for β/γ ≲ 1, individuals with connectivity ≳10-fold above average are likely to seed epidemics. These rare events correspond to infections seeded in a connectivity hub in the network. Conversely, individuals that have a connectivity smaller than the average connectivity in the population have an extinction probability close to 1 in the same regime. In the next section, we will use Q_p(R₀, k₀) as a weighting factor for testing protocols.

Zoom In Zoom Out Reset image size

A similar pattern is observed in the conditional proliferation probability Q_p(R₀, n), which depends on an intermediate outbreak size n as defined above (figure 4(a)). This similarity is intuitively clear: an individual of degree k₀ generates, after one transmission step, an infected cluster of size n proportional to k₀.

In the remainder of this section, we sketch the derivation of the above results from the outbreak statistics on networks (more details are given in methods; readers interested primarily in the application to contact protocols may skip this part). First, following reference [13], we define the basic building block of the epidemic dynamics on a network: the probability T that an infected individual transmits the infection to a given contact; in other words, the probability that the infection spreads along a given branch of the network. This probability determines the reproductive rate on the network,

$$\begin{equation}{R}_{0}(T)=\langle \!\langle (k-1)\,T\,\rangle \!\rangle =\frac{\left\langle k(k-1)\right\rangle }{\left\langle k\right\rangle }\,T.\end{equation} \tag{ 10 }$$

Here ⟨⋯⟩ denotes averaging over the degree distribution $\hat{p}(k)$ and 《⋯》 averaging over the nearest-neighbor degree distribution ${\hat{p}}_{2}(k)=k\hat{p}(k)/\left\langle k\right\rangle$. The resulting growth rate in the exponential regime, λ(R₀, γ, σ) is then given by equation (5).

The branch transmission probability T, which has been used as an independent parameter in previous work, has no direct analogue in a well-mixed system. However, we can express T in terms of the rate parameters β and γ,

$$\begin{equation}T(\beta ,\gamma )=\frac{1}{1+\frac{\langle k\rangle }{\beta /\gamma }},\end{equation} \tag{ 11 }$$

where we have used equation (9) (see methods). Together with equation (10), we obtain

$$\begin{equation}{R}_{0}(\beta ,\gamma )=\left(\frac{\left\langle {k}^{2}\right\rangle -\left\langle k\right\rangle }{\langle {k\rangle }^{2}+\frac{\beta }{\gamma }\left\langle k\right\rangle }\right)\frac{\beta }{\gamma };\end{equation} \tag{ 12 }$$

this relation compares the reproductive number on a contact network with the corresponding well-mixed system, which has R₀ = β/γ. It displays the two network effects discussed above: preferential spreading on high-degree nodes and saturation of susceptible contacts (figure 3(b)). In the well-mixed model limit, (⟨k⟩ ≫ 1 with Var(k) = o(⟨k⟩²)), we get R₀ → β/γ as expected.

Next, we use the branch transmission probability T to compute the probability that an outbreak proliferates, given that it originates on a node of degree k₀. Again, a previous result for the SIR model [13, 20] carries over to the SEIR model, because incubation delays growth but does not affect the probability of transmission. The proliferation probability takes the form

$$\begin{equation}{Q}_{\text{p}}({R}_{0},{k}_{0})=1-{\left[1-T({R}_{0})+{q}_{b}T({R}_{0})\right]}^{{k}_{0}}.\end{equation} \tag{ 13 }$$

Here q_b is the probability that a sub-outbreak originating on a given branch of the network does not proliferate; this probability can be obtained from a self-consistent summation procedure (reference [13], see methods). The product on the rhs of equation (13) says that transmissions along each of the k₀ branches originating from patient zero are statistically independent sub-outbreaks. The relation between T and R₀ is given by equation (10).

The global proliferation probabilities are then given by averages over the contact network [13, 20]. We obtain

$$\begin{equation}{Q}_{\text{p}}({R}_{0})=\left\langle {Q}_{\text{p}}({R}_{0},k)\right\rangle \end{equation} \tag{ 14 }$$

and

$$\begin{equation}{Q}_{\text{p}}({R}_{0},n)=1-{\left\langle \!\left\langle 1-{Q}_{\text{p}}({R}_{0},k-1)\right\rangle \!\right\rangle }^{n}.\end{equation} \tag{ 15 }$$

The product on the rhs of equation (15) says that transmissions originating from each of the n initially infected individuals are statistically independent sub-outbreaks; this form is analogous to equation (13). A given individual has been infected by one of its contacts and has (k − 1) contacts left to transmit the infection further. We note that these expressions take into account the contacts of typical individuals in the outbreak rather than individual contacts of specific individuals; for this reason, the average is taken using the nearest-neighbor degree distribution ${\hat{p}}_{2}(k)$ [13]. As before, the conditional proliferation probability determines the establishment size of epidemics in a contact network, n*(R₀), by the condition Q_p(R₀, n*) = 1/2. Below, we use this function to characterize the performance of surveillance testing protocols.

The structures of the epidemiological model and of the contact network have important implications for the analysis of empirical data, in particular, for predicting the likelihood of future local outbreaks.

First, the reproductive number R₀ is often inferred from the observed growth rate λ. In this procedure, we have to take into account infection and incubation periods, as given by the SEIR relation in equation (5). Using the corresponding SIR relation, which neglects incubation, the reproductive number can be drastically underestimated (figure 5(a)).

Zoom In Zoom Out Reset image size

Second, given a correctly inferred value of R₀, the contact network shapes the corresponding proliferation probability of local outbreaks, as given by equations (13)–(15). Using the naive expression for a well-mixed system, Q_p = 1/R₀, proliferation can be substantially overestimated (figure 5(b)); the same applies to the conditional probabilities Q_p(R₀, k₀) and Q_p(R₀, n). The correct inference of R₀ and of proliferation probabilities is a crucial step in the surveillance and containment protocols discussed below.

As mentioned above, the degree of the seeding node has a big impact on the probability of having an epidemic. Using Bayes' rule, we calculate the probability $\hat{p}({k}_{0}\vert \mathrm{p}\mathrm{r}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n})$ that an epidemic originates from a node of degree k₀,

$$\begin{equation}\hat{p}({k}_{0}\vert \mathrm{p}\mathrm{r}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n})\sim \frac{{Q}_{\text{p}}({R}_{0},{k}_{0})}{{Q}_{\text{p}}({R}_{0})}.\end{equation} \tag{ 16 }$$

We upgraded the testing surveillance by using the weight function $\omega ({k}_{0})\equiv \hat{p}({k}_{0}\vert \mathrm{p}\mathrm{r}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n})$ for the random sampling of the protocol. With this weighting, we test highly connected individuals more often than poorly connected ones (figure A2 in methods).

We find an appreciable improvement of the degree-based testing protocol compared to the uniform testing protocol. The performance is better both in terms of the cluster size of the outbreak (figure 7) and in the time of detection (see figure A3 in methods). Consistently, the improvement is smaller for larger values of R₀, where the effects of the network are less significant.

Zoom In Zoom Out Reset image size

Figure 7 shows the ratio between the cluster size at the moment when the first infected individual is detected in the degree-based sampling protocol, n_k , and in the uniform sampling protocol, n_U. Figure 7(b) shows that when the population is well-mixed, there is no significant difference between the two protocols. However, when the network of contacts regulates the interactions between individuals, as shown in figure 7(a), the expected cluster size is reduced in all the range of parameters that we explore. This reduction is greater for smaller values of R₀, reaching values larger than 20%. We also measure a greater improvement for smaller sample sizes in the protocol, consistent with the fact that with larger uniform random samples it is more likely to detect infectious hubs that drive the proliferation.

We use a classical epidemiological model consisting of susceptible (S), exposed (E), infected (I) and recovered (R) individuals. Exposed individuals are infected individuals that do not test positive and cannot infect other individuals yet. The total population size is N. Here we consider the deterministic dynamics of the system described by the ODE system of equations (1)–(4). We are interested in the first stage of the dynamics, where the approximation S ≈ N is valid. With this approximation, and assuming the initial conditions S(t) = N − 1, E(t) = R(t) = 0 and I(t) = 1, equations (2) and (3) can be written as

$$\begin{equation}\dot{E}(t)=-\sigma E+\beta I\end{equation} \tag{ A.1 }$$

$$\begin{equation}\dot{I}(t)=\sigma E-\gamma I.\end{equation} \tag{ A.2 }$$

The solution of the system of equation (A.2) is given by

$$\begin{equation}\left(\begin{array}{c}E(t)\\ I(t)\end{array}\right)=C{e}^{{\lambda }_{+}t}\left(\begin{array}{c}\quad 1\\ \frac{{\lambda }_{+}+\sigma }{\beta }\end{array}\right)-C{e}^{{\lambda }_{-}t}\left(\begin{array}{c}\quad 1\\ \frac{{\lambda }_{-}+\sigma }{\beta }\end{array}\right),\end{equation} \tag{ A.3 }$$

where the constant C is determined by the initial condition and λ₊ and λ₋ are the eigenvalues of the matrix

$$\begin{equation}\left(\begin{matrix}\hfill -\sigma \hfill & \hfill \beta \hfill \\ \hfill \sigma \hfill & \hfill -\gamma \hfill \end{matrix}\right)\end{equation} \tag{ A.4 }$$

that are given by

$$\begin{equation}{\lambda }_{\pm }=-\frac{\sigma +\gamma }{2}\pm \frac{1}{2}\sqrt{{(\sigma -\gamma )}^{2}+4\sigma \beta }.\end{equation} \tag{ A.5 }$$

As soon as the outbreak is seeded, the time evolution of E and I is lead by the first term in equation (A.3) because λ₋ < 0. With the initial conditions E(0) = 0 and I(0) = 1, equation (A.3) can be written as

$$\begin{equation}\left(\begin{matrix}\hfill E(t)\hfill \\ \hfill I(t)\hfill \end{matrix}\right)=\frac{\beta }{\sqrt{{(\sigma -\gamma )}^{2}+4\sigma \beta }}\left(\begin{matrix}\hfill 1\hfill \\ \hfill \frac{{\lambda }_{+}+\sigma }{\beta }\hfill \end{matrix}\right){e}^{{\lambda }_{+}t}.\end{equation} \tag{ A.6 }$$

As mentioned above, we decided to include the exposed compartment of the population because it affects dramatically the dynamics of the infected individuals (figure 1(b)). More concretely, for SARS-Cov-2 the estimated average time that takes an exposed individual to become infected τ_σ is in the same order of magnitude of other relevant time-scale, namely the average infectious period τ_γ [1]. Additionally, under these circumstances the classical relation between the exponential growth rate and the parameter R₀, the initial reproductive number, defined as the quotient between the infectious rate and the recovery rate β/γ, need to be revisited.

The approximation S ≈ N is only valid when few individuals are infected and breaks down when the infection approaches its expected peak. From equations (1)–(4), it is possible to derive an expression for the number of infected individuals in the peak. Using the condition $\dot{I}=0$ and the constrain of constant total population size N = S + E + I + R, we have that

$$\begin{equation}{I}_{\text{peak}}\approx N\cdot \frac{\sigma \gamma }{\beta (\sigma +\gamma )}\left(\beta /\gamma -1-\mathrm{log}\,\beta /\gamma \right).\end{equation} \tag{ A.7 }$$

We expect then deterministic exponential grow when I(t) < I_peak. In figure 2(a), for the case of β/γ = 4.0, the peak is reached at I_peak ≈ 1210 individuals, and for that reason we see that the average epidemic trajectory grows exponentially for the whole range shown. On the other hand, in figure 2(b) for the case β/γ = 1.3, the peak is reached at I_peak ≈ 86. Therefore, an exponential growth is observed only for a shorter period and the approximation made before is invalid for the later part.

We generate the random network of contacts by following a preferential attachment method [12]. Such a method produces graphs with heavy-tailed degree distribution $\hat{p}(k)\sim {k}^{-\alpha }$, like the one shown in figure A1(a). Under these circumstances, most of the individuals have few number of contacts close to the mean degree and a non-negligible amount of individuals have a large number of contacts. The degree of the highly connected individuals can reach values up to two orders of magnitude larger than the mean degree. In such case, a finite variance is guaranteed by a cutoff in the degree distribution at large degree values as a consequence of the finite size of the network.

Zoom In Zoom Out Reset image size

In this context, a new variable T_nm has been defined in the literature [13] to characterize the epidemiological dynamics on the random network, which accounts for the probability that the infection is transmitted through an edge in the network between nodes n and m. The average value, T, is given in terms of the typical effective time of duration of an infection τ, which is a exponentially distributed stochastic variable with mean value τ_γ = γ⁻¹, and the rate of infectious encounters between individuals in the network ρ. Given that in our model ρ is proportional to β (see equation (9)), we derived that for the case of a SEIR dynamics, T is given by

$$\begin{align}T& =1-{\int }_{0}^{\mathrm{\infty }}\gamma {\mathit{\text{e}}}^{-\gamma {\tau }^{\mathrm{\prime }}}{e}^{-\frac{\beta }{\langle k\rangle }{\tau }^{\mathrm{\prime }}}\;\mathrm{d}{\tau }^{\mathrm{\prime }}\\ & =1-\frac{1}{1+\frac{\beta /\gamma }{\langle k\rangle }},\end{align} \tag{ A.8 }$$

where the second term correspond to the average probability across the network that no such transmission event occurs.

Moreover, a critical value ${T}_{\text{c}}=\left\langle k\right\rangle /(\left\langle {k}^{2}\right\rangle -\left\langle k\right\rangle )$ has been derived for an uncorrelated network previously, which sets the minimum value of T required for an outbreak to have the chance to turn into an epidemic. The value of T_c corresponds to the correction factor in equation (10) that we derived for the renormalized basic reproductive number in the network of contacts.

As it has been derived before [13], we start the derivation with a Dyson-like self-consistent equation for the probability q_b that an individual at the end of a randomly selected interacting edge will not produce an epidemic,

$$\begin{equation}{q}_{b}=\frac{{\sum }_{k=1}^{\infty }k\,\hat{p\,}(k){\Psi}(T,k-1)}{\left\langle k\right\rangle },\end{equation} \tag{ A.9 }$$

where

$$\begin{equation}{\Psi}(T,k)={\left(1-T+(T{q}_{b})\right)}^{k}\end{equation} \tag{ A.10 }$$

and 1 − T is the probability that the disease does not spread through the interacting edge, Tq_b is the probability that it spreads through the edge but it does not produce an epidemic and ${\hat{p}}_{2}(k)=k\hat{p}(k)/\left\langle k\right\rangle$ is the distribution of interacting edges connected to an individual at the end of a randomly selected interacting edge. This nearest-neighbors degree distribution is important because by following a randomly chosen interacting edge, is more likely to reach highly connected hubs than by randomly selecting single nodes. We solved equation (A.9) numerically for a fixed value of T. Here, we are able to write the proliferation probability equations in terms of the epidemiological parameters that define the dynamics of equations (1)–(4) and (6)–(8) for the case of the SEIR model or by the epidemiological measurable values λ, τ_γ and τ_σ .

Near the critical value R₀ = 1, some properties of the system show well defined scaling behavior [21, 41]. In particular, the proliferation probability ${Q}_{\text{p}}\sim {(T-{T}_{\text{c}})}^{-{\beta }_{\text{p}}}$ with β_c = 1/(3 − α) when 2 < α < 3, which is the case for the networks generated by the preferential attachment method and that corresponds to the most common networks occurring in nature [41]. In our case, some approximations may become invalid due to large fluctuations. Given the relation that we derive between T and R₀, it might be possible to derive new scaling laws for R₀ approaching its corresponding critical value. However, this is beyond the scope of this work.

Zoom In Zoom Out Reset image size

The function ω(k) defined in the main text is proportional to the likelihood function L(k) = p(proliferation|k) = Q_p(R₀, k). It should be reduced to a constant value when the network of contacts is turned off because then Q_p(R₀, k) = Q_p(R₀)

Figure A3(a) shows the outcomes of the simulations for the average time of first detection to complement the evaluation of the performance of the degree-based sampling testing protocol. Similar to the cluster size shown in figure 7(a), there is a significant improvement in the time at which the first infected individual is detected. For small values of R₀ the reduction of such time can be large than 20%.

Zoom In Zoom Out Reset image size

To test that the random testing protocol has an effect on a structured population characterized by a network of contacts, we simulate the degree-based sampling protocol turning off the network of interaction as a negative control but still using the likelihood function L described above for the sampling of individuals. Figures 7(b) and A3(b) show the performance of the control protocol in terms of cluster size and time of the first detection, respectively. In both cases the performance of the degree-based protocol is the same as the uniform sampling protocol.

The simulations were performed with a modified version of the individual-based SEIR python code seirsplus publicly available in https://github.com/ryansmcgee/seirsplus. The simulations consist of a Gillespie algorithm for independent Poisson processes.