Intervention Threshold for Epidemic Control in Susceptible-Infected-Recovered Metapopulation Models

TL;DR

This paper analyzes how to determine the critical intervention threshold in a metapopulation SIR model to prevent global epidemics, emphasizing targeted strategies over random interventions.

Contribution

It provides a theoretical framework for the intervention threshold in SIR metapopulation models and demonstrates the effectiveness of targeted interventions through simulations.

Findings

Targeted intervention to high-degree patches is more effective.

Critical fraction of low-risk patches prevents global outbreaks.

Theoretical results are validated by Monte Carlo simulations.

Abstract

Metapopulation epidemic models describe epidemic dynamics in networks of spatially distant patches connected with pathways for migration of individuals. In the present study, we deal with a susceptible-infected-recovered (SIR) metapopulation model where the epidemic process in each patch is represented by an SIR model and the mobility of individuals is assumed to be a homogeneous diffusion. Our study focuses on two types of patches including high-risk and low-risk ones, in order to evaluate intervention strategies for epidemic control. We theoretically analyze the intervention threshold, indicating the critical fraction of low-risk patches for preventing a global epidemic outbreak. We show that targeted intervention to high-degree patches is more effective for epidemic control than random intervention. The theoretical results are validated by Monte Carlo simulation for synthetic and realistic scale-free patch networks. Our approach is useful for exploring better local interventions aimed at containment of epidemics.

Equations

$$d_{k{k}^{\prime }}$$

$$N_k$$

$${\lambda }_{k{k}^{\prime }}$$

$$D_k^n=\sum _{k^{\prime }}{D}_{k^{\prime }}^{n-1}(k^{\prime }-1)P(k|k^{\prime })\left(1-\frac{D_k^{n-1}}{V_k}\right)\left(1-R_0^{-{\lambda }_{k^{\prime }k}}\right),$$

$${\lambda }_{k^{\prime }k}$$

$$P(k|k^{\prime })$$

$$1-R_0^{-{\lambda }_{k^{\prime }k}}$$

$$D_k^n$$

$${\Theta }^n=\frac{p\alpha \bar{N}}{\mu }\frac{⟨k^2⟩-⟨k⟩}{⟨k{⟩}^2}(R_0-1){\Theta }^{n-1}.$$

$$R_{\ast }$$

$$D_{k,H}^n$$

$$D_{k,L}^n$$

$${\lambda }_{k^{\prime }k}^H$$

$${\lambda }_{k^{\prime }k}^L$$

$$1-(R_0^H{)}^{-{\lambda }_{k^{\prime }k}}$$

$$1-(R_0^L{)}^{-{\lambda }_{k^{\prime }k}}$$

$$\left(1-\frac{D_{k,H}^{n-1}}{V_{k,H}}\right)\simeq 1\mathrm{and}\left(1-\frac{D_{k,L}^{n-1}}{V_{k,L}}\right)\simeq 1.$$

$${\Theta }_H^n$$

$${\Theta }_L^n$$

$$[k^{\alpha }]$$

$${\Theta }_H^n$$

$${\Theta }_L^n$$

$$\displaystyle\left(\begin{array}[]{c}\Theta_{H}^{n}\\ \Theta_{L}^{n}\end{array}\right)=J\left(\begin{array}[]{c}\Theta_{H}^{n-1}\\ \Theta_{L}^{n-1}\end{array}\right),$$

$$\displaystyle J:=\left(\begin{array}[]{cc}\frac{p\bar{N}\alpha_{H}(R_{0}^{H}-1)}{\mu}(\phi_{1}-\phi_{2})&\frac{p\bar{N}\alpha_{L}(R_{0}^{H}-1)}{\mu}(\phi_{1}-\phi_{2})\\ \frac{p\bar{N}\alpha_{H}(R_{0}^{L}-1)}{\mu}\phi_{2}&\frac{p\bar{N}\alpha_{L}(R_{0}^{L}-1)}{\mu}\phi_{2}\end{array}\right),$$

$${\varphi }_1:=(⟨k^2⟩-⟨k⟩)/⟨k{⟩}^2,$$

$${\varphi }_2:=([k^2]-[k])/⟨k{⟩}^2.$$

$$\frac{p\bar{N}}{\mu }\{{\alpha }_H(R_0^H-1)({\varphi }_1-{\varphi }_2)+{\alpha }_L(R_0^L-1){\varphi }_2\}.$$

$$R_c$$

$$[k^{\alpha }]$$

$$0\le q(k)\le 1,$$

$${\int }_{k}q(k)p(k)dk$$

$$q^{\mathrm{rn}}(k)$$

Full text

††preprint: APS/123-QED

I Introduction

In the modern age of expanding globalization, epidemic spreading is a serious matter of global public health. Countermeasures, such as vaccinations, antiviral medications, and social distancing, have been practiced for controlling past infectious diseases. However, emerging and re-emerging infectious diseases pose perpetual challenges for controlling them due to environmental changes and diversification of human behavior Fauci and Morens (2012); Morens and Fauci (2013). Therefore, it is significant to continuously explore systematic methods for planning effective epidemic control strategies. Mathematical models are powerful tools for understanding epidemic spreading processes which are complex phenomena involved in the type of disease, host immunity, environmental conditions, and human mobility patterns Keeling and Rohani (2011). Mathematical methods have been widely used to estimate epidemic outcomes and evaluate the effectiveness of preventive measures.

There are a variety of mathematical models for epidemic spreading, from simple to complex ones. Compartment epidemic models assuming homogeneous mixing of individuals are classical and simple Diekmann and Heesterbeek (2000). These models have been extended to more complex and realistic ones by incorporating additional factors, such as social structures Strang (1991), spatial structures Riley et al. (2015), seasonal forcing Tanaka and Aihara (2013); Bjørnstad and Viboud (2016); Buonomo et al. (2018), human mobility patterns Barrat et al. (2008); Balcan et al. (2010); Belik et al. (2011a); Bajardi et al. (2011); Tizzoni et al. (2014); Urabe et al. (2016). Metapopulation epidemic models are a class of models that describe epidemic spreading processes in a group of spatially separated patches connected via migration pathways. This framework has been intensively studied to develop analytical methods for global epidemic thresholds Sattenspiel and Dietz (1995); Keeling et al. (2004); Colizza and Vespignani (2007, 2008); Belik et al. (2011b); Balcan and Vespignani (2011); Wang et al. (2012); Apolloni et al. (2014); Wang et al. (2014); Tizzoni et al. (2015); Wang et al. (2017); Gómez-Gardeñes et al. (2018); Granell and Mucha (2018); Soriano-Paños et al. (2018) and widely applied to understanding actual epidemic outbreaks Cador et al. (2016); De Luca et al. (2018); Lee et al. (2018); Meakin et al. (2019); Laager et al. (2019). In a metapopulation model, infection and recovery events occur in each patch and migration of individuals potentially causes global epidemic spreading. Metapopulation models have often been employed to consider inhomogeneous mixing of individuals. The SIR metapopulation model with fully connected patches was analyzed to examine the properties of the global basic reproduction number (which is differentiated from the local basic reproduction number in an isolated patch) governing the global epidemic threshold Arino and Van den Driessche (2003); Allen and Van den Driessche (2008). The global reproduction number can be numerically estimated using the next-generation matrix approach, but a theoretical issue is to derive its explicit expression as a function of system parameters including local epidemiological parameters and human mobility patterns Cross et al. (2005).

Colizza and Vespignani Colizza and Vespignani (2007, 2008) derived an analytical expression of the global invasion threshold (i.e. the global reproduction number) for an SIR metapopulation model with complex patch connectivity under several assumptions. They clarified the effect of heterogeneous network connectivity on the global epidemic threshold. The analysis was conducted under the assumption that the local reproduction number is the same for all the patches. However, the conditions of patches are thought to be heterogeneous in reality. In fact, it was reported that the local reproduction numbers estimated from real data for seasonal influenza are different between local areas Chowell et al. (2010). The heterogeneity of local reproduction numbers can be partially attributed to the difference in the immunization coverage rates in local areas. Under the patch heterogeneity, the effectiveness of strategic interventions for epidemic control has been evaluated using a susceptible-infected-susceptible (SIS) metapopulation model in our previous study Tanaka et al. (2014). The result shows that targeted interventions for high-degree patches are more effective than random interventions. However, it is still unclear whether this result holds for other types of metapopulation models. For instance, the SIR epidemic process representing an epidemic outbreak as a transient state is qualitatively different from the SIS epidemic process representing an endemic state as a stationary state. In fact, theoretical approaches are completely different between them; the global epidemic threshold of an SIR metapopulation model is analyzed based on a branching process method Colizza and Vespignani (2007, 2008), while that of an SIS metapopulation model is analyzed based on local stability of the disease-free equilibrium state Colizza et al. (2007); Tanaka et al. (2014). Thus, the impact of local intervention in SIR metapopulation models with patch heterogeneity still remains to be studied.

In the present study, we aim to analyze the intervention threshold in SIR metapopulation models consisting of high-risk and low-risk patches as shown in Fig. 1. This model framework is similar to that considered in Ref. Tanaka et al. (2014), but the epidemic process is qualitatively different. We introduce an intervention rate $u$ representing the fraction of low-risk patches that have received intervention. When $u=0$, all the patches are high-risk and a global epidemic outbreak inevitably occurs. In the other extreme case with $u=1$, all the patches are low-risk and a global epidemic outbreak is prevented. Therefore, we can expect that there is a certain critical value $u=u_{c}\in(0,1)$ (called an intervention threshold), separating the outbreak and non-outbreak regimes. We use this threshold as a measure to compare different intervention strategies. The smaller the intervention threshold is, the more effective the intervention strategy is. The main novelty of this study is to theoretically derive the intervention threshold $u_{c}$ for random and targeted interventions. Furthermore, our theoretical results are validated by numerical simulations. The comparison of the intervention thresholds shows that targeted intervention for high-degree patches is more effective than random intervention. Our result indicating the effectiveness of targeted intervention in metapopulations reminds us of the effectiveness of targeted immunization in complex contact networks of individuals Pastor-Satorras and Vespignani (2002); Pastor-Satorras et al. (2015). However, these model frameworks are largely different because dynamical processes within network nodes are considered in metapopulation models but not in contact network models.

In Sec. II, we first introduce the analysis framework proposed in the previous study Colizza and Vespignani (2007, 2008) and then describe our approach. In Sec. III, we show theoretical and numerical results. In Sec. IV, we conclude this study.

II Methods

II.1 SIR metapopulation model with identical patches

We first introduce a method for analyzing a global epidemic threshold in an SIR metapopuation model with identical patches, following Refs. Colizza and Vespignani (2007, 2008). We will extend this method to the case with non-identical patches in the subsequent section.

An SIR metapopulation model describes epidemic spreading in a network of spatially separated patches, interconnected with migration pathways. The number of patches is denoted by $V$. The patches are assumed to totally contain a sufficiently large number $N$ of individuals, who are susceptible (S), infected (I), or recovered (R). Epidemic dynamics in each patch follows an SIR process Kermack and McKendrick (1932): a susceptible individual changes to an infected one with transmission rate $\beta$ when contacting with an infected individual ($\mbox{S}+\mbox{I}\rightarrow 2\mbox{I}$); an infected individual changes to a recovered one with recovery rate $\mu$ ($\mbox{I}\rightarrow\mbox{R}$). We assume homogeneous mixing of individuals in each patch where the local reproduction number is given by $R_{0}=\beta/\mu$. To allow global epidemic spreading, the local reproduction number $R_{0}$ needs to be larger than unity. Individuals can migrate from one patch to a neighboring one through the pathway. This is regarded as a diffusion process Colizza et al. (2007) and the diffusion rate from a patch is denoted by $p$. Under a homogeneous diffusion process, the diffusion rate from a patch with degree $k$ to one of the neighboring patches (with any degree $k^{\prime}$) is given by

[TABLE]

In a stationary state, the number of individuals in a patch with degree $k$ is obtained as follows Colizza and Vespignani (2007, 2008):

[TABLE]

where $\bar{N}=N/V$ is the average population size per patch and $\langle k\rangle$ is the mean degree.

We consider an initial condition that an infected individual invades a metapopulation system of susceptible individuals. The total number of infected individuals in a patch is proportional to the stationary population in the patch, described as $\alpha N_{k}$ for a patch with degree $k$, where the coefficient $\alpha$ depends on the type of disease and other factors. The average infection period of an infected individual is given by the inverse of the recovery rate, $\mu^{-1}$. Therefore, if an epidemic occurs in a patch with degree $k$, the average number of infected individuals who move to a neighboring patch with degree $k^{\prime}$ is represented as follows Colizza and Vespignani (2007, 2008):

[TABLE]

We focus on the time evolution of the number of “infected” patches which are defined as the patches that undergo an outbreak. The analysis is based on the basic branching process Ball et al. (1997); Harris (2002). We denote by $D_{k}^{0}$ the number of infected patches with degree $k$ at generation 0 (i.e. in the beginning of the process). These patches bring about new infected patches with degree $k$ in their neighborhood, the number of which is represented as $D^{1}_{k}$ at generation 1. In this way, we define as $D^{n}_{k}$ the number of infected patches with degree $k$ at generation $n$. Assuming that the number of infected patches is sufficiently small in the early stage of the process, we can approximately relate $D^{n}_{k}$ to $D^{n-1}_{k}$ as follows Colizza and Vespignani (2007, 2008):

[TABLE]

where $P(k)$ denotes the degree distribution of the patch network and $V_{k}$ denotes the number of patches with degree $k$. This equation is derived based on the notion that each infected patch with degree $k^{\prime}$ at generation $(n-1)$ will spread infection in the $(k^{\prime}-1)$ neighboring patches except the one that originally transmitted infection, the probability that a neighboring patch of a patch with degree $k^{\prime}$ has degree $k$ is $P(k|k^{\prime})$, and the probability that the disease does not become extinct when $\lambda_{k^{\prime}k}$ infected individuals invade in a patch with $R_{0}$ is given by $(1-R_{0}^{-\lambda_{k^{\prime}k}})$ Bailey et al. (1975); Murray (2005).

For analytical tractability, we deal with special cases under the assumptions of homogeneous diffusion in mobility, uncorrelated patch networks (i.e. without degree-degree correlation), and the local reproduction number close to an epidemic threshold. From Eqs. (1)-(3), the number of seeds of infection is given by

[TABLE]

In an uncorrelated patch network, the following equation holds Dorogovtsev and Mendes (2002):

[TABLE]

When the local reproduction number is close to the epidemic threshold, i.e. $R_{0}-1\ll 1$, the outbreak probability is approximated as follows:

[TABLE]

By substituting Eqs. (5)-(7) into Eq. (II.1), we obtain the following equation Colizza and Vespignani (2007, 2008):

[TABLE]

By defining $\Theta^{n}:=\sum_{k^{\prime}}D^{n}_{k^{\prime}}(k^{\prime}-1)$, the above equation can be rewritten by the following recurrence formula:

[TABLE]

The condition that $\Theta^{n}$ does not increase with $n$ is given by Colizza and Vespignani (2007, 2008)

[TABLE]

where $R_{*}$ represents the global reproduction number. If $R_{0}$ is close to 1, then $\alpha\simeq 2(R_{0}-1)/(R_{0})^{2}$ according to Ref. Murray (2005). Using this approximation, Eq. (10) is simplified as follows:

[TABLE]

II.2 SIR metapopulation model with high-risk and low-risk patches

Extending the framework in Sec. II.1, we analyze an SIR metapopulation model consisting of high-risk and low-risk patches for examining epidemic intervention strategies Allen et al. (2007); Tanaka et al. (2014). We assume that only a fraction of patches can receive intervention and become low-risk due to budgetary constraints. The local reproduction number in such low-risk patches is denoted by $R_{0}^{L}$ and that in the remaining high-risk patches is by $R_{0}^{H}(>R_{0}^{L})$.

Let us define $D^{n}_{k,H}$ and $D^{n}_{k,L}$ as the numbers of infected high-risk and low-risk patches with degree $k$ at generation $n$, respectively. The numbers of individuals who experience the disease during an outbreak in the high-risk and low-risk patches are represented as $\alpha_{H}N_{k}$ and $\alpha_{L}N_{k}$, respectively. The numbers of seeds from high-risk and low-risk patches with degree $k$ are denoted by $\lambda_{k^{\prime}k}^{H}$ and $\lambda_{k^{\prime}k}^{L}$, respectively. We define $Q(k)$ as the probability that a randomly chosen patch with degree $k$ is a low-risk one. Considering the transmission of infection from high-risk and low-risk patches separately, the recurrence formulae for $D_{k,H}^{n}$ and $D_{k,L}^{n}$ are written as follows (as in Eq. (4)):

[TABLE]

[TABLE]

where $V_{k,H}$ and $V_{k,L}$ represent the numbers of high-risk and low-risk patches with degree $k$, respectively.

As in Eq. (5), we obtain

[TABLE]

Assuming $R_{0}^{H}\simeq 1$ and $R_{0}^{L}\simeq 1$, we can use the following approximations (as in Eq. (7)):

[TABLE]

In the early stage of the propagation, it follows

[TABLE]

By defining

[TABLE]

we can rewrite Eqs. (12)-(13) as follows:

[TABLE]

These recurrence equations are simply written as follows:

[TABLE]

where

[TABLE]

The eigenvalues of $J$ are given by 0 and

[TABLE]

The condition that $\Theta^{n}_{H}$ and $\Theta^{n}_{L}$ do not diverge in the limit of $n\to\infty$ is equivalent to the condition that the absolute values of all the eigenvalues of $J$ are smaller than 1. Hence, the condition that a global outbreak does not occur is given by:

[TABLE]

where $R_{c}$ represents the global reproduction number in the case that high-risk and low-risk patches coexist. When $R_{0}$ is close to unity, $\alpha=2(R_{0}-1)/R_{0}^{2}$. Therefore, the global reproduction number is rewritten as follows:

[TABLE]

where $\psi(x):=2(x-1)^{2}/x^{2}$. Based on this formula, we can evaluate the intervention threshold for different strategies as described in Sec. III.1.

II.3 Numerical simulation methods

We describe numerical methods for simulating epidemic propagation processes in SIR metapopulation models, which are used to validate our theoretical results. The state of each individual is susceptible (S), infected (I), or recovered (R). Initially the population in patch $j$ is set at $N_{j}=k_{j}\bar{N}/\langle k\rangle$ for $j=1,\ldots,V$. The numbers of susceptible, infected, and recovered individuals in patch $j$ are denoted by $S_{j}$, $I_{j}$, and $R_{j}$, respectively. A patch whose degree is close to $\langle k\rangle$ is chosen to have ten initial infected individuals. The remaining individuals are susceptible. We consider discrete-time dynamical processes and denote the unitary time step by $\tau$. At each time step, the state of each individual in patch $j$ is probabilistically updated. The update process consists of two stages: epidemic and mobility stages. In the epidemic stage, each susceptible individual turns into an infected one with probability $1-(1-\beta_{j}\tau/N_{j})^{I_{j}}$ and each infected individual turns into a recovered one with probability $\mu\tau$. After all individuals have been updated in the epidemic stage, the mobility stage starts. In the mobility stage, each individual moves to one of the neighboring patches with probability $p\tau$. The above procedure is repeated for all the individuals at each time step and continued for finite time steps until infected individuals disappear.

We set $\bar{N}=1000$, $\tau=0.1$, and $\mu=1$ for all the patches, $\beta_{j}=2$ for high-risk patches, and $\beta_{j}=1.01$ for low-risk patches, unless otherwise noted. To focus on epidemic spreading in heterogeneous patch networks, we employed synthetic scale-free networks with $V=200$ patches having degree distribution $P(k)\sim k^{-\gamma}$ with $\gamma=2.1$ generated by the uncorrelated configuration model Catanzaro et al. (2005) and the real US airport network having a scale-free property, containing $V=500$ patches Colizza et al. (2007). We performed 50 simulations with different network realizations for each parameter condition.

III Results

First, we theoretically analyze the intervention threshold in SIR metapopulation models with heterogeneously connected patches in Sec. III.1. We deal with random and targeted interventions Tanaka et al. (2014). Then, we numerically validate the theoretical results in Sec. III.2.

III.1 Theoretical results

For theoretical analysis, we approximate the degree $k$ as a continuous variable by representing its expectation value over many realizations of networks Barabási et al. (2016). We consider a probability density function $p(k)$ for a continuous degree distribution, instead of the discrete degree distribution $P(k)$. We also define a probability density function $q(k)$ for a continuous intervention probability, instead of the discrete intervention probability $Q(k)$. Accordingly, the summations with respect to $k$ in the previous section are replaced with integrals over $k$. In particular, we redefine the brackets in Eq. (21) as follows:

[TABLE]

The total intervention rate $u$ represents the fraction of low-risk patches. For a given $u$, we need to appropriately define $q(k)$ such that

[TABLE]

III.1.1 Threshold for random intervention

First, we deal with random intervention, where the low-risk patches are chosen at random. Namely, the probability that a patch is low-risk is constant independently of the patch degree. From Eq. (40), we obtain the probability density function $q(k)$ for the random intervention as follows:

[TABLE]

In this case, we have $[k]=u\langle k\rangle$ from Eq. (38) and $\phi_{2}=u\phi_{1}$ from Eqs. (33), (34), and (38). Using these equations and Eq. (36), the global reproduction number $R_{c}$ is described as follows:

[TABLE]

By solving $R_{c}^{\rm rn}=1$ with respect to $u$, we obtain the critical intervention threshold as follows:

[TABLE]

above which a global epidemic outbreak is prevented.

The global reproduction number $R_{c}^{\rm rn}$ for random intervention in a scale-free patch network is computed from Eq. (LABEL:eq:R_random) and plotted as a function of the intervention rate $u$ and the mobility rate $p$ in Fig. 2. The yellow filled circles represent the critical intervention threshold $u^{\rm rn}_{c}$ given by Eq. (43). As seen from Fig. 2, $R_{c}^{\rm rn}$ decreases monotonically with the intervention rate $u$ and increases monotonically with the mobility rate $p$. The fact that the value of $u^{\rm rn}_{c}$ increases with $p$ suggests that more interventions for epidemic control are required when spatial movements of individuals are more active.

The theoretical result in Eq. (43) also reveals the influence of the underlying mobility network on the intervention threshold. The connectivity of the scale-free patch network is varied with the degree exponent $\gamma$, which is typically in the range $2<\gamma<3$ Barabási et al. (2016). For a scale-free network with a degree distribution $p(k)\sim k^{-\gamma}$, the maximum degree $k_{\rm max}$ and the minimum one $k_{\rm min}$ satisfy the following relationship Barabási et al. (2016):

[TABLE]

Therefore, a smaller value of $\gamma$ means a larger difference between the maximum and minimum degrees. The intervention threshold given by Eq. (43) is plotted as a function of the degree exponent $\gamma$ and the mobility rate $p$ in Fig. 3. We can see that, for a fixed value of the mobility rate $p$, the intervention threshold $u_{c}$ increases with a decrease in $\gamma$. This result shows that more interventions are required for preventing a global outbreak in scale-free patch networks having a hub patch with a larger maximum degree. At $\gamma\simeq 2$, we have $k_{\rm max}\simeq V$ from Eq. (44) and there exists a hub patch that connects to almost all other patches. Such a big hub patch easily causes a global outbreak and leads to a large critical intervention rate.

III.1.2 Threshold for targeted intervention

Next, we consider targeted intervention, which means that important patches are preferentially selected to be low-risk. Here we measure the importance of a patch using the degree centrality Freeman (1978); the more connections a patch has, the more likely it is chosen as a low-risk patch. In this case, $q(k)$ should be a monotonically increasing function of $k$.

As a candidate of such a function, we define a piecewise function $q_{l}(k)$ as shown in Fig. 4(a), represented as follows:

[TABLE]

where $l$ is a real value ranging between $k_{\rm min}$ and $k_{\rm max}$. We define the expectation value of $q_{l}(k)$ with respect to $k$ as follows:

[TABLE]

For scale-free networks with $p(k)\sim k^{-\gamma}$, we can show that $\hat{q}_{l}$ is monotonically decreasing with increasing $l$. In the limit of $l\rightarrow k_{\rm min}$, $\hat{q}_{l}$ approaches the maximum value 1. When $l=k_{\rm max}$, $\hat{q}_{l}$ takes the minimum value $\hat{q}_{k_{\rm max}}$.

We define the probability density function $q(k)$ separately for the two cases of $u\geq\hat{q}_{k_{\rm max}}$ and $u<\hat{q}_{k_{\rm max}}$. If $u\geq\hat{q}_{k_{\rm max}}$, we can find $k^{*}\in[k_{\rm min},k_{\rm max}]$ such that $\hat{q}_{k^{*}}=u$, satisfying Eq. (40). Therefore, we use $q_{k^{*}}$ as $q(k)$. Otherwise, the piecewise function with any $l$ does not satisfy Eq. (40). In this case, we use a non-piecewise function as shown in Fig. 4(b). The probability density function $q(k)$ for the targeted intervention is defined as follows:

[TABLE]

We can show that the latter case also satisfies the requirements for the probability density function, Eqs. (39)-(40), as follows:

[TABLE]

[TABLE]

For the probability density function $q(k)=q^{\rm tg}(k)$, the global reproduction number $R^{\rm tg}_{c}$ is obtained as Eq. (37) which depends on $\phi_{2}$ in Eq. (34) where $[k^{2}]$ and $[k]$ are computed from Eq. (38) with Eq. (47). The global reproduction number $R_{c}^{\rm tg}$ for a scale-free patch network is shown as a function of the intervention rate $u$ and the mobility rate $p$ in Fig. 5. As in the random intervention case, $R_{c}^{\rm tg}$ increases with decreasing $u$ and increasing $p$. We can obtain the critical intervention threshold $u^{\rm tg}_{c}$ by numerically solving $R_{c}^{\rm tg}=1$ with respect to $u$. The yellow filled circles represent the values of $u^{\rm tg}_{c}$ which increase with $p$. By comparing Fig. 5 with Fig. 2, it can be visually confirmed that the critical intervention threshold for targeted intervention ($u^{\rm tg}_{c}$) is much smaller than that for random intervention ($u^{\rm rn}_{c}$).

The critical intervention threshold $u^{\rm tg}_{c}$ is obtained by solving $R_{c}=1$ with respect to $u$ using Eq. (37) and Eq. (47). The threshold is plotted as a function of the degree exponent $\gamma$ and the mobility rate $p$ in Fig. 6. In contrast to the random intervention case, there is no remarkable change in $u_{c}$ with increasing $\gamma$ for a fixed value of $p$. This result implies that a very small intervention rate with which a few hub patches are changed to low-risk ones is sufficient for preventing global outbreaks.

III.1.3 Comparison of the intervention thresholds

The global reproduction number $R_{c}$ in Eq. (37) is different between the random and targeted interventions, because $\phi_{2}$ depends on $q(k)$. A smaller value of $R_{c}$ for the same intervention rate $u$ means a more effective intervention strategy. We show that the targeted intervention is more effective than the random one. It is sufficient to prove the following inequality:

[TABLE]

where $\phi_{2}^{\rm rn}$ and $\phi_{2}^{\rm tg}$ denote $\phi_{2}$ (Eq. (34)) for $q(k)=q^{\rm rn}(k)$ and $q(k)=q^{\rm tg}(k)$, respectively. We first deal with the case of $u\geq\hat{q}_{k_{\rm max}}$ and then that of $u<\hat{q}_{k_{\rm max}}$.

First, we assume $u\geq\hat{q}_{k_{\rm max}}$. We can evaluate $\Delta\phi_{2}(u)$ as follows:

[TABLE]

From Eq. (45), the last term is equivalent to

[TABLE]

Now we define the following functions:

[TABLE]

where $b(k)$ is a monotonically increasing function of $k$ for a scale-free network, satisfying $b(k_{\rm min})<0$, $b(k_{\rm max})>0$, and $\int_{k}b(k)dk=0$. Using these functions, we can represent $\Delta\phi_{2}(u)$ as follows:

[TABLE]

Therefore, the inequality (48) holds if $B_{k^{*}}$ is non-negative. From Eq. (50), we have

[TABLE]

From the monotonicity of $b(k)$, $dB_{l}/dl$ in Eq. (52) is a unimodal function. Therefore, we have

[TABLE]

Hence, $B_{l}$ is a monotonically increasing function of $l$. From $B_{k_{\rm min}}=0$, we obtain $B_{k^{*}}\geq 0$. From Eq. (51), $\Delta\phi_{2}(u)\geq 0$ is satisfied.

Next, we assume $u<\hat{q}_{k_{\rm max}}$. From Eq. (48), it follows

[TABLE]

From $k\geq k_{\rm min}>0$, we have $k(k-\langle k\rangle)^{2}\geq 0$, yielding $\langle k(k-\langle k\rangle)^{2}\rangle\geq 0$. This yields $\langle k^{3}\rangle\geq 2\langle k\rangle\langle k^{2}\rangle-\langle k\rangle^{3}\geq 0$. From this property, we can evaluate Eq. (55) as follows:

[TABLE]

Therefore, $\Delta\phi_{2}(u)\geq 0$ holds. From Eq. (55), we find that $\Delta\phi_{2}(u)$ is a monotonically increasing function of $u$.

III.2 Numerical validation

We numerically study the effect of the local interventions on the final epidemic size. The final epidemic size is measured by the ratio of individuals who have experienced the disease during an outbreak period, given by $R_{\infty}/N$ where $R_{\infty}$ equals to $\sum_{j}R_{j}$ after the outbreak. Due to the finiteness of the number of degrees in simulations, we used a discretized version of Eq. (47) for the targeted intervention.

The average values of the final epidemic size over 50 simulations are plotted against the intervention rate $u$ in Fig. 7(a) for synthetic scale-free patch networks generated with the configuration model Catanzaro et al. (2005) and in Fig. 7(b) for the U.S. airport network representing the connectivity of flight routes between 500 major airports in United States Colizza et al. (2007). We see that, in both networks, the targeted intervention is much more effective than the random intervention as it requires a much lower intervention rate for containment of epidemics.

In Fig. 8, the numerical results of the final epidemic size are shown for variation of the intervention rate $u$ and the mobility rate $p$. Figures 8(a) and (b) correspond to the results for random and targeted interventions in synthetic scale-free patch networks, respectively. A comparison between these two figures obviously shows that the targeted intervention is more effective than the random intervention for reducing the epidemic size. The same property is confirmed in Figs. 8(c) and (d), which correspond to random and targeted interventions in the US airport network, respectively. In all the cases, the final epidemic size decreases with increasing $u$ and increases with increasing $p$. The theoretical values of the intervention threshold $u_{c}$ are superimposed as yellow filled circles, indicating that they are in good agreement with the thresholds which are recognized from the numerical results.

IV Conclusion and Discussion

In the present study, we have analyzed the intervention threshold in SIR metapopulation models with scale-free patch connectivity, consisting of high-risk and low-risk patches. Under the assumption that a high-risk patch is changed to a low-risk patch by reducing the local reproduction number by an intervention, we have compared the effectiveness of random and targeted interventions through theoretical and numerical analyses. The theoretical results have shown that the intervention targeted to high-degree patches is more effective than the random intervention. They have been validated by the numerical simulations using the synthetic scale-free networks and the realistic US airport network. Our result indicating the effectiveness of the targeted intervention for SIR metapopulation models is consistent with that for SIS metapopulation models in a similar framework Tanaka et al. (2014). As the global reproduction number is expressed as a function of the intervention rate and the mobility rate, one can calculate the critical intervention threshold for a given mobility rate and estimate the minimum effort for containment of epidemics. We have found that a higher human mobility rate leads to a larger intervention threshold. This finding suggests that travel restrictions are effective, especially when using targeted intervention. We have also revealed that the intervention threshold is larger for a scale-free patch network with a smaller degree exponent. This result implies that an existence of a very big hub patch increases the difficulty of preventing global outbreaks.

The framework for examining intervention strategies in this study has a potential to be extended to more realistic cases in terms of human mobility patterns and intervention strategies. There are other types of human mobility patterns, such as recurrent (commuting) mobility Belik et al. (2011b); Balcan and Vespignani (2011); Panigutti et al. (2017); Gómez-Gardeñes et al. (2018); Granell and Mucha (2018) and adaptive mobility Meloni et al. (2011); Wang et al. (2012). Moreover, human mobility networks can be better estimated from higher-resolution data such as real-world traffic network data Merler and Ajelli (2009); Balcan et al. (2010), mobile phone data Wesolowski et al. (2012); Tizzoni et al. (2014); Panigutti et al. (2017), and GPS data Vazquez-Prokopec et al. (2013). It is also intriguing to test other intervention strategies, such as those based on other network centralities, because the important patches are not necessarily high-degree ones. Another strategy is to combine the intervention to local patches and travel restrictions. It would be possible that the optimal intervention strategy is different depending on how to evaluate the epidemic outcome. Therefore, an appropriate assessment of the social impact of global epidemics from microscopic and macroscopic levels is becoming increasingly important Ball et al. (2015); Massaro et al. (2018).