Optimal Allocation of Resources for Suppressing Epidemic Spreading on Networks

TL;DR

This paper develops an optimal resource allocation strategy to suppress epidemic spreading on networks, showing that proportional allocation to node degree maximizes epidemic threshold and that in strong infection regions, low-degree nodes require more resources.

Contribution

It introduces a theoretical framework for optimal resource distribution in epidemic control, revealing counterintuitive strategies in strong infection scenarios.

Findings

Optimal curing rate proportional to node degree maximizes epidemic threshold.

In strong infection regions, low-degree nodes should receive more resources.

Theoretical results are validated through simulated annealing.

Abstract

Efficient allocation of limited medical resources is crucial for controlling epidemic spreading on networks. Based on the susceptible-infected-susceptible model, we solve an optimization problem as how best to allocate the limited resources so as to minimize the prevalence, providing that the curing rate of each node is positively correlated to its medical resource. By quenched mean-field theory and heterogeneous mean-field (HMF) theory, we prove that epidemic outbreak will be suppressed to the greatest extent if the curing rate of each node is directly proportional to its degree, under which the effective infection rate $\lambda$ has a maximal threshold $\lambda_c^{opt}=1/\left\langle k \right\rangle$ where $\left\langle k \right\rangle$ is average degree of the underlying network. For weak infection region ($\lambda\gtrsim\lambda_c^{opt}$), we combine a perturbation theory with Lagrange multiplier method (LMM) to derive the analytical expression of optimal allocation of the curing rates and the corresponding minimized prevalence. For general infection region ($\lambda>\lambda_c^{opt}$), the high-dimensional optimization problem is converted into numerically solving low-dimensional nonlinear equations by the HMF theory and LMM. Counterintuitively, in the strong infection region the low-degree nodes should be allocated more medical resources than the high-degree nodes to minimize the prevalence. Finally, we use simulated annealing to validate the theoretical results.

Equations

$$⟨{\mu }_i⟩=\mu and{\mu }_i⩾0,\forall i.$$

$$\frac{d{\rho }_i(t)}{dt}=-{\mu }_i{\rho }_i(t)+\beta \left[1-{\rho }_i(t)\right]\sum _j{A}_{ij}{\rho }_j(t).$$

$${\rho }_i=\frac{\beta {\sum }_j{A}_{ij}{\rho }_j}{{\mu }_i+\beta {\sum }_j{A}_{ij}{\rho }_j}.$$

$${\Lambda }_{max}\ge {\text{v}}^T(-\mathbb{U}+\beta \mathbb{A})\text{v},$$

$${\Lambda }_{max}\ge {\text{v}}^T(-\bar{\mathbb{U}}+\beta \mathbb{A})\text{v}-{\text{v}}^T\Delta \mathbb{U}\text{v}=-\mu +\beta ⟨k⟩.$$

$${\mu }_i={\mu }_i^{\ast }=\mu \frac{k_i}{⟨k⟩},$$

$${\lambda }_{c,QMF}^{opt}=\frac{1}{⟨k⟩}.$$

$$⟨{\mu }_k⟩={\sum }_{k}P(k){\mu }_k=\mu and{\mu }_k⩾0,\forall k,$$

$$\frac{d{\rho }_k\left(t\right)}{dt}=-{\mu }_k{\rho }_k\left(t\right)+\beta \left[1-{\rho }_k\left(t\right)\right]k\Theta (t),$$

$$\Theta (t)=\sum _k\frac{kP(k)}{⟨k⟩}{\rho }_k(t).$$

$${\rho }_k=\frac{\beta k\Theta }{{\mu }_k+\beta k\Theta }.$$

$$\Theta =\sum _k\frac{kP(k)}{⟨k⟩}\frac{\beta k\Theta }{{\mu }_k+\beta k\Theta }.$$

$${\beta }_{c,HMF}=\frac{⟨k⟩}{{\sum }_k\frac{k^{2}P(k)}{{\mu }_k}}.$$

$$\mathcal{L}=\sum _k\frac{k^{2}P(k)}{{\mu }_k}+\tau \left(\sum _{k}P(k){\mu }_k-\mu \right),$$

$$\frac{\partial \mathcal{L}}{\partial {\mu }_k}=-\frac{k^{2}P(k)}{{\mu }_k^2}+\tau P(k),$$

$${\lambda }_{c,HMF}^{opt}=\frac{1}{⟨k⟩},$$

$${\mu }_k={\mu }_k^{\ast }=\mu \frac{k}{⟨k⟩}$$

$$\Delta \Theta =\frac{1}{{\beta }^2⟨k⟩}\sum _k\frac{P(k)}{k}\Delta {\mu }_k^2.$$

$$\Delta \rho =-\frac{{\Theta }^{\ast }}{\beta }\sum _k\frac{P(k)}{k}\Delta {\mu }_k+\left(1-{\Theta }^{\ast }\right)\Delta \Theta .$$

$$\Delta {\rho }^{opt}$$

$$\Delta {\mu }_k=\frac{\mu }{2}⟨k⟩\lambda \left(⟨k⟩-k\right)\Delta \lambda \simeq \frac{\mu }{2}\left(⟨k⟩-k\right)\Delta \lambda .$$

$${\mu }_k=\{\begin{array}{c}\sqrt{\frac{\beta k\Theta }{\tau }+\frac{\kappa \beta k^2}{\tau ⟨k⟩}}-\beta k\Theta >0,k<k_c\text{hfill}\\ 0,k\ge k_c\text{hfill}\end{array}$$

$$\sqrt{\frac{\beta \tau ⟨k⟩}{\Theta }}\sum _{k=k_{\min }}^{k_{\max }}\sqrt{\xi }P(k)-\beta \tau ⟨k⟩\sum _{k=k_{\min }}^{k_{\max }}\xi P(k)$$

$$-\beta \kappa \tau \sum _{k=k_{\min }}^{k_{\max }}k\xi P(k)-\frac{\kappa }{⟨k⟩\Theta }\sum _{k=k_c}^{k_{\max }}kP(k)=0,$$

$$\mu =\sqrt{\frac{\beta \Theta }{\tau ⟨k⟩}}\sum _{k=k_{\min }}^{k_c}k{\xi }^{-\frac{1}{2}}P(k)-\beta \Theta \sum _{k=k_{min}}^{k_c}kP(k),$$

$$\Theta =\sqrt{\frac{\beta \tau \Theta }{⟨k⟩}}\sum _{k=k_{\min }}^{k_c}{k}^2{\xi }^{-\frac{1}{2}}P(k)+\frac{1}{⟨k⟩}\sum _{k=k_c}^{k_{\max }}kP(k),$$

$$\frac{\beta }{⟨k⟩}\sum _k\frac{k^{2}P(k)}{{\mu }_k+\beta k\Theta }=1.$$

$$\begin{array}{rl}& \sum _k{\frac{\partial f}{\partial {\mu }_k}|}_{\left({\mu }_k^{\ast },{\Theta }^{\ast }\right)}\Delta {\mu }_k+{\frac{\partial f}{\partial \Theta }|}_{\left({\mu }_k^{\ast },{\Theta }^{\ast }\right)}\Delta \Theta \\ & +\frac{1}{2}\sum _k{\sum _{k^{\prime }}\frac{{\partial }^{2}f}{\partial {\mu }_k\partial {\mu }_{k^{\prime }}}|}_{\left({\mu }_k^{\ast },{\Theta }^{\ast }\right)}\Delta {\mu }_k\Delta {\mu }_{k^{\prime }}\\ & +\sum _k{\frac{{\partial }^{2}f}{\partial {\mu }_k\partial \Theta }|}_{\left({\mu }_k^{\ast },{\Theta }^{\ast }\right)}\Delta {\mu }_k\Delta \Theta +{\frac{1}{2}\frac{{\partial }^{2}f}{\partial \Theta \partial \Theta }|}_{\left({\mu }_k^{\ast },{\Theta }^{\ast }\right)}\Delta {\Theta }^2=0,\end{array}$$

$$\begin{array}{rl}& {\frac{\partial f}{\partial {\mu }_k}|}_{\left({\mu }_k^{\ast },{\Theta }^{\ast }\right)}=-\frac{P(k)}{\beta ⟨k⟩}\\ & {\frac{\partial f}{\partial \Theta }|}_{\left({\mu }_k^{\ast },{\Theta }^{\ast }\right)}=-1\\ & {\frac{{\partial }^{2}f}{\partial {\mu }_k\partial {\mu }_{k^{\prime }}}|}_{\left({\mu }_k^{\ast },{\Theta }^{\ast }\right)}={\delta }_{k{k}^{\prime }}\frac{2P(k)}{{\beta }^2⟨k⟩k}\\ & {\frac{{\partial }^{2}f}{\partial {\mu }_k\partial \Theta }|}_{\left({\mu }_k^{\ast },{\Theta }^{\ast }\right)}=\frac{2P(k)}{\beta ⟨k⟩}\\ & {\frac{{\partial }^{2}f}{\partial \Theta \partial \Theta }|}_{\left({\mu }_k^{\ast },{\Theta }^{\ast }\right)}=2.\end{array}$$

$$\begin{array}{rl}& -\frac{1}{\beta ⟨k⟩}\sum _{k}P(k)\Delta {\mu }_k-\Delta \Theta +\frac{1}{{\beta }^2⟨k⟩}\sum _k\frac{P(k)}{k}\Delta {\mu }_k^2\\ & +\frac{2}{\beta ⟨k⟩}\sum _{k}P(k)\Delta {\mu }_k\Delta \Theta +\Delta {\Theta }^2=0.\end{array}$$

$$\Delta \Theta =\frac{1}{{\beta }^2⟨k⟩}\sum _k\frac{P(k)}{k}\Delta {\mu }_k^2.$$

Full text

A challenging problem in epidemiology is how best to allocate limited resources of treatment and vaccination so that they will be most effective in suppressing or reducing outbreaks of epidemics. This problem has been a subject of intense research in statistical physics and many other disciplines Pastor-Satorras et al. (2015); Nowzari et al. (2016). Inspired by the percolation theory, the simplest strategy is to randomly choose a fraction of nodes to immunize. However, the random immunization is inefficient for heterogeneous networks. Later on, many more effective immunization strategies have been developed, ranging from global strategies like targeted immunization based on node degree Pastor-Satorras and Vespignani (2002) or betweenness centrality Holme et al. (2002) to local strategies, like acquaintance immunization Cohen et al. (2003) and (bias) random walk immunization Holme (2004); Stauffer and Barbosa (2006) and to some others in between Gomez-Gardenes et al. (2006). Further improvements were done by graph partitioning Chen et al. (2008) and the optimization of the susceptible size Schneider et al. (2011). Besides the degree heterogeneity, community structure has also a major impact on disease immunity Masuda (2009); Salathé and Jones (2010). Recently, a message-passing approach was used to find an optimal set of nodes for immunization Altarelli et al. (2014). The immunization has been mapped onto the optimal percolation problem Morone and Makse (2015). Based on the idea of explosive percolation, an “explosive immunization” method has been proposed Clusella et al. (2016). However, some diseases like the common cold and influenza that can be modeled by the susceptible-infected-susceptible (SIS) model, do not confer immunity and individuals can be infected over and over again. Under the situations, one way to control the spread of the diseases is to reduce the risk of the infection, such as adaptive rewiring links incident to infected individuals Gross et al. (2006) and dynamical interplay between awareness and epidemic spreading Granell et al. (2013).

An alternative way to control the epidemic spreading of SIS type by designing an optimal strategy for distributing the limited medical resources so as to suppress the epidemic outbreak to the greatest extent and minimize the prevalence once the epidemic outbreak has happened. It is reasonable to assume the curing rate of each node is positively correlated to the medical resources allocated to it. Therefore, the optimal allocation of medical resources is equivalent to that of the curing rates. Assuming the total medical resources are limited, the average curing rate is thus considered to be fixed. This problem has been addressed as a constraint optimization problem in several previous works. When the curing rate can be only tuned in a fixed number of feasible values, this problem has been proved to be NP-complete Prakash et al. (2013). Instead, when the curing rate can continuously varies in a given interval, some efficient algorithms have been developed for minimizing the threshold of epidemic outbreak Wan et al. (2008); Preciado et al. (2013) or the steady-state infection density Gourdin et al. (2011). In the present work, we theoretically solve the constraint optimization problem in both epidemic-free and endemic phases within the mean-field framework. On the one hand, we prove that the epidemic outbreak can be suppressed to the most extent when the curing rate of each node is directly proportional to its degree, under which the epidemic threshold is maximized that is the inverse of the average degree of the underlying network. On the other hand, once the epidemic has broken out but close to the threshold, we analytically show the optimal curing rate should be adjusted in terms of the difference of node degree with average degree and the distance to epidemic threshold. For the general infection region, the optimization problem can be simplified to solve three nonlinear equations.

To formulate our problem, we consider the SIS model on an undirected network of size $N$. The network is described by an adjacency matrix $\mathbb{A}$ whose entries are defined as $A_{ij}=1$ if nodes $i$ and $j$ are connected, and $A_{ij}=0$ otherwise. Each node is either susceptible or infected. A susceptible node $i$ can be infected by its infective neighbor with an infection rate $\beta$, and an infected node $i$ recovers with a nonvanishing curing rate $\mu_{i}$. Here, we consider that the curing rate is allowed to vary from one node to another one. In general, the more available medical resource of a node $i$ has, the larger $\mu_{i}$ is. Assuming that the total amount of medicine resource is limited, the average curing rate is thus fixed, i.e.,

[TABLE]

Our goal is to find out an optimal allocation of $\{\mu_{i}\}$ under the constraint Eq.(1) so as to minimize the prevalence $\rho$, that is the fraction of infected nodes.

In the quenched mean-field (QMF) theory, the probability $\rho_{i}(t)$ that node $i$ is infected at time $t$ is described by $N$-intertwined equations Y. Wang ; Mieghem et al. (2009); Gómez et al. (2010),

[TABLE]

In the steady state, $d\rho_{i}(t)/dt=0$, $\rho_{i}$ is determined by a set of nonlinear equations,

[TABLE]

One can notice that $\rho_{i}=0$ is always a solution of Eq.(2). This trivial solution corresponds to an absorbing state with no infective nodes. A nonzero solution $\rho_{i}>0$ exists if the effective infection rate $\lambda=\beta/\mu$ is larger than the so-called epidemic threshold $\lambda_{c}$. In this case, the prevalence $\rho=\sum\nolimits_{i}{{\rho_{i}}}/N$ is nonzero corresponding to an endemic state. By linear stability analysis for Eq.(2) around $\rho_{i}=0$, $\lambda_{c}$ is determined by which the largest eigenvalue of the matrix, $-\mathbb{U}+\beta\mathbb{A}$, is zero, where $\mathbb{U}=diag(\mu_{i})$ is a diagonal matrix. For the standard SIS model, $\mu_{i}\equiv\mu$ for all $i$, one can immediately obtain the well-known result, $\lambda_{c,QMF}^{sta}=1/\Lambda_{max}(\mathbb{A})$ with the largest eigenvalue of the adjacency matrix $\Lambda_{max}(\mathbb{A})$. In our SIS model, the outbreak of epidemics will be suppressed to the greatest extent, which implies that the epidemic threshold of the optimal SIS model will be maximized.

For this purpose, we first decompose the diagonal matrix $\mathbb{U}$ into two diagonal matrices, $\mathbb{U}=\mathbb{\bar{U}}+\Delta\mathbb{U}$, where $\mathbb{\bar{U}}=diag\{\mu k_{i}/\left\langle k\right\rangle\}$ with $k_{i}$ being the degree of node $i$ and $\Delta\mathbb{U}=diag\{\Delta\mu_{i}\}$. Since $Tr(\mathbb{U})=Tr(\mathbb{\bar{U}})=N\mu$, $\Delta\mathbb{U}$ must satisfy the constraint $Tr(\Delta\mathbb{U})=0$. For the real symmetric matrix, $\mathbb{U}-\beta\mathbb{A}$, its largest eigenvalue $\Lambda_{max}$ satisfies the following inequality,

[TABLE]

where v is a column vector satisfying $\textbf{v}\in\mathbb{R}^{N}$ and $||\textbf{v}||=1$. If we set $\textbf{v}=\frac{1}{{\sqrt{N}}}{\left({1,\cdots,1}\right)^{T}}$, Eq.(4) becomes

[TABLE]

Since $\Lambda_{max}=0$ at the epidemic threshold, Eq.(5) leads to an upper bound of epidemic threshold, $\lambda_{c}\leq 1/\left\langle k\right\rangle$. The condition that the epidemic threshold equals to the upper bound holds when v is the eigenvector of $\mathbb{U}-\beta\mathbb{A}$ corresponding to its largest eigenvalue. If we set $\mathbb{U}=\mathbb{\bar{U}}$ and $\beta=\mu/\left\langle k\right\rangle$, $-\mathbb{U}+\beta\mathbb{A}=-\mu/\left\langle k\right\rangle\mathbb{L}$, where $\mathbb{L}$ is the Laplacian matrix of the underlying network. It is well-known that the smallest eigenvalue of $\mathbb{L}$ is zero and the corresponding eigenvector is v. Therefore, if the curing rate of each node is directly proportional to its degree, i.e.,

[TABLE]

the epidemic threshold will be maximized,

[TABLE]

In the QMF theory, the epidemic threshold of the optimal SIS model is no less than that of the standard SIS model, $\lambda_{c,QMF}^{opt}\geq\lambda_{c,QMF}^{sta}$, as the lower bound of $\Lambda_{max}(\mathbb{A})$ is $\left\langle k\right\rangle$ for any types of networks Mieghem (2011).

The above results can be also derived from the heterogeneous mean-field (HMF) theory. In the framework of HMF, these nodes with the same degree are considered to be statistically equivalent. The constraint Eq.(1) becomes

[TABLE]

where $\mu_{k}$ is the curing rate of nodes of degree $k$, and $P(k)$ is the degree distribution. Some related works have studied the SIS model Dezső and Barabási (2002) and its metapopulation version Shen et al. (2012) with the curing rate, $\mu_{k}\sim k^{\alpha}$, but such a power-law form did not guarantee to be the optimal one. In Borgs et al. (2010), the authors consider a simple heuristic strategy to control epidemic extinction where the curing rate is directly proportional to node degree. They showed that on any graph with bounded degree the extinction time is sublinear with the size of the network. Further improvement has been done by a heuristic PageRank algorithm to allocate curing rates based on the initial condition of infected nodes Chung et al. (2009). The present study does not require any assumptions about the form of the cure rate with node degree in advance except to the constraint Eq.(8). The dynamical evolution of $\rho_{k}(t)$, the probability of nodes of degree $k$ being infected at time $t$, reads Pastor-Satorras and Vespignani (2001),

[TABLE]

where $\Theta$ is the probability of finding an infected node following a randomly chosen edge. In the case of uncorrelated networks, $\Theta(t)$ can be written as

[TABLE]

In the steady state, $d\rho_{k}(t)/dt=0$, Eq.(9) becomes

[TABLE]

Substituting Eq.(11) into Eq.(10), we obtain a self-consistent equation of $\Theta$,

[TABLE]

The epidemic threshold is determined by which the derivation of the r.h.s of Eq.(12) with respect to $\Theta$ at $\Theta=0$ equals to one, leading to

[TABLE]

For a given $P(k)$, maximizing $\beta_{c}$ is equivalent to minimizing the denominator of the r.h.s of Eq.(13). For this purpose, we employ Lagrange multiplier method (LMM) to maximize the epidemic threshold, where the Lagrange function is written as,

[TABLE]

where $\tau$ is called the Lagrange multiplier. Taking the derivation of $\mathcal{L}$ with respect to $\mu_{k}$,

[TABLE]

and letting ${{\partial\mathcal{L}}\mathord{\left/{\vphantom{{\partial L}{\partial{\mu_{k}}=0}}}\right.\kern-1.2pt}{\partial{\mu_{k}}=0}}$ combined with Eq.(8), we arrive at a maximal epidemic threshold

[TABLE]

and the corresponding allocation of $\{\mu_{k}\}$,

[TABLE]

Interestingly, the HMF results are consistent with the QMF ones. Also, in the HMF theory the epidemic threshold of the optimal SIS model is no less than that of the standard SIS model, $\lambda_{c,HMF}^{opt}\geq\lambda_{c,HMF}^{sta}=\left\langle k\right\rangle/\left\langle k^{2}\right\rangle$.

For $\lambda$ larger than but close to $\lambda_{c}^{opt}$, $\lambda\gtrsim\lambda_{c}^{opt}$, we shall combine a perturbation theory with LMM to optimize the prevalence. To the end, we assume that for $\lambda=\lambda_{c}^{opt}+\Delta\lambda$, $\mu_{k}=\mu_{k}^{*}+\Delta\mu_{k}$ and $\Theta=\Theta^{*}+\Delta\Theta$, where $\Theta^{*}=1-\frac{\mu}{{\beta\left\langle k\right\rangle}}$ is the solution of Eq.(12) for $\mu_{k}=\mu_{k}^{*}$. Expanding Eq.(12) around ${\left({\mu_{k}^{*},{\Theta^{*}}}\right)}$ to the second order, and then using the constraint $\sum\nolimits_{k}P(k)\Delta\mu_{k}=0$ and simultaneously ignoring the second-order small quantity $\Delta\Theta^{2}$, it yields [See Appendix A for details]

[TABLE]

Around ${\left({\mu_{k}^{*},{\Theta^{*}}}\right)}$, the change $\Delta\rho$ in the prevalence $\rho=\sum\nolimits_{k}P(k)\rho_{k}$ can be written as

[TABLE]

Again using LMM to minimize $\Delta\rho$ under the constraints $\sum\nolimits_{k}P(k)\Delta\mu_{k}=0$ and Eq.(18), we obtain a minimal $\rho=\rho^{*}+\Delta\rho^{opt}$ with ${\rho^{*}}=\sum\nolimits_{k}{P(k)\frac{{\beta k{\Theta^{*}}}}{{\mu_{k}^{*}+\beta k{\Theta^{*}}}}}$ and

[TABLE]

Since $\left\langle{{k^{-1}}}\right\rangle>{\left\langle k\right\rangle^{-1}}$ for any degree inhomogeneous networks in terms of Jensen’s inequality, $\Delta\rho^{opt}<0$ and thus $\rho$ will be reduced. The corresponding optimal allocation $\mu_{k}=\mu_{k}^{*}+\Delta\mu_{k}$ with

[TABLE]

This implies that as $\lambda$ is increased from $\lambda_{c}^{opt}$, the curing rates of the nodes with degrees less than the average degree will be increased, while the curing rates of the nodes with degrees larger than the average degree will be decreased. The amplitude of the change will depend on the difference between the degree of each node and the average degree, $\left\langle k\right\rangle-k$, and the distance of the effective infection rate to its critical value, $\Delta\lambda$.

For $\lambda$ is larger than but not close to $\lambda_{c}^{opt}$, $\lambda>\lambda_{c}^{opt}$, since the nonlinear characteristic of the model, analytical expression of optimal allocation of $\{\mu_{k}\}$ and the corresponding the minimal $\rho$ is almost impossible. However, with the aid of HMF theory and LMM, the high-dimensional optimization problem can be converted to numerically solving the low-dimensional nonlinear equations [See Appendix B for details]. In the general infection region, $\mu_{k}$ satisfies the following equation,

[TABLE]

where $\tau$ and $\kappa$ are the Lagrange multipliers, and $k_{c}$ is a threshold degree to guarantee $\mu_{k}>0$ for $k<k_{c}$ and it will be determined later. $\Theta$, $\tau$ and $\kappa$ are determined by the following three equations,

[TABLE]

[TABLE]

[TABLE]

where we have used $\xi=k/({\left\langle k\right\rangle+\kappa k})$.

To numerically solve $\Theta$, $\tau$ and $\kappa$ by Eqs.(26,27,28), $k_{c}$ is needed to be known in advance. To the end, we adopt a numerical scheme as follows. (i) Firstly we set $k_{c}=k_{max}$ where $k_{max}$ is the maximal degree of the underlying network; (ii) we numerically solve $\Theta$, $\tau$ and $\kappa$ by Eqs.(26,27,28), and then test the condition $\mu_{k}>0$ for all $k<k_{c}$ by Eq.(25); (iii) if the condition is not satisfied, $k_{c}$ will be decreased by $k_{c}\leftarrow k_{c}-1$ and return to ii) until the condition Eq.(25) is fulfilled.

Figure 1 shows the optimized results of $\rho$ as a function $\lambda$ (solid line) in Erdös-Rényi (ER) random networks (a) and Barabási-Albert scale-free networks (b) with equal network size $N=1000$ and average degree $\left\langle k\right\rangle=4$. For comparison, we also show the results of the standard SIS model (dotted line) and of the SIS model with the curing rates $\mu_{i}=\mu_{i}^{*}$ (dashed line). As expected by the theoretical prediction, the epidemic threshold of the optimal SIS model $\lambda_{c}^{opt}=1/\left\langle k\right\rangle$, which is significantly larger than that of the standard SIS model, but coincides with the case of $\mu_{i}=\mu_{i}^{*}$. While for $\lambda>\lambda_{c}^{opt}$, the prevalence for $\mu_{i}=\mu_{i}^{*}$ is always larger than the optimal choice, and even larger than the standard SIS model in the strong infection region, indicating that $\mu_{i}=\mu_{i}^{*}$ is not a good choice once the epidemic outbreak has happened.

We use the simulated annealing (SA) technique to validate our theoretical results. The SA builds a Monte Carlo Markov Chain that in the long run converges to the minimum of a given energy function $\mathcal{E}$, where $\mathcal{E}=\rho$ can be obtained by numerically iterating Eq.(3). The main steps of SA are as follows. At beginning, we assign to a given set of $\{\mu_{i}\}$ satisfying the constraint Eq.(1) (e.g., $\mu_{i}=\mu$ for all $i$). Then, we randomly choose two distinct nodes, say $i$ and $j$, and try to make the changes $\mu_{i}\leftarrow\mu_{i}+\delta$ and $\mu_{j}\leftarrow\mu_{j}-\delta$ with the standard Metropolis probability $\min(1,e^{-\beta_{SA}\Delta\mathcal{E}})$, where $\delta$ is randomly chosen between $-\mu_{i}$ and $\mu_{i}+\mu_{j}$ to guarantee the curing rate is always not less than zero. $\beta_{SA}$ is the inverse temperature of SA which slowly increases from $10^{-2}$ to $10^{4}$ via an annealing protocol. $\Delta\mathcal{E}$ is the change of the energy function $\mathcal{E}$ due the change of $\mu_{i}$ and $\mu_{j}$, We tested several different annealing protocols and we adopted one in which the inverse temperature of SA $\beta_{SA}$ is updated by $\beta_{SA}\leftarrow 1.01\beta_{SA}$ after each $N$ attempts for updating $\{\mu_{i}\}$. The SA results are also shown in Fig.1 (square dots), which agree with the theoretical prediction.

In Fig.2, we show the optimal allocation of $\{\mu_{k}\}$ as a function of node degree $k$ for several distinct $\lambda$ in ER random networks (a) and BA scale-free networks (b), in which the theoretical results and the SA ones are indicated by the lines and dots, respectively. For $\lambda\gtrsim\lambda_{c}^{opt}$, $\mu_{k}$ increases linearly as $k$ with the slope depending on the distance to the epidemic threshold. The results have been well predicted by Eq.(17) and Eq.(21). For the region away from the threshold, $\mu_{k}$ will deviate from linear relation with $k$. For sufficiently large $\lambda$, $\mu_{k}$ for large $k$ can be less than that for small $k$, and even $\mu_{k}$ vanishes when $k$ exceeds a threshold value, as given by Eq.(25). This surprising result implies that in the strong infection region more medicine resources should be put into these low-degree nodes other than high-degree nodes.

In conclusion, we have theoretically studied a constraint optimization problem as how best to distribute the limited medicine resources (curing rates) for controlling the epidemics of SIS type. Based on the QMF and HMF theories, we have shown that the optimal allocation lies in the effective infection rate $\lambda$ (or the basic reproduction number $R_{0}=\left\langle k\right\rangle\lambda$). If $R_{0}\leqslant 1$, the curing rate of each node should be in direct proportion to its degree, under which the epidemic outbreak will be suppressed to the most extent and the epidemic threshold will be maximized, Eq.(7) or Eq.(16). Once the maximal epidemic threshold is just across ($R_{0}\gtrsim 1$), the epidemic will spread persistently. In this case, we have analytically shown that the change in the curing rate of each node depends linearly on the difference between the average degree and its degree and the distance to epidemic threshold, Eq.(21). For the general infection region ($R_{0}>1$), it is almost impossible to derive an analytical solution of the optimization problem; however, it can be simplified to an much more easily problem of numerical calculation of three nonlinear equations, Eqs.(26,27,28). Surprisingly, we found that in the strong infection region the curing rates of the low-degree nodes can overpower those of the high-degree nodes to ensure the minimization of the prevalence.

An interesting generalization is how to solve the present constraint optimization problem based on other existing theoretical methods, such as pair mean-field method that takes into account the role of dynamical correlations between neighboring nodes Eames and Keeling (2002); Gleeson (2011); Boguñá et al. (2013); Mata et al. (2014); Kiss et al. (2015); Cai et al. (2016); Cator and Van Mieghem (2012); Mata and Ferreira (2013); I. Z. Kiss (2015). Moreover, the method presented here could be applied to a number of other optimization problems, for example, controlling opinion dynamics in social networks Castellano et al. (2009). This will be the subject of future work.

Acknowledgments: This work was supported by National Science Foundation of China (Grants Nos. 11205002, 61473001, 21673212), the Key Scientific Research Fund of Anhui Provincial Education Department (Grant No. KJ2016A015) and “211” Project of Anhui University (Grant No. J01005106).

Appendix A Weak infection region

For $\lambda$ larger than but close to $\lambda_{c}^{opt}$, $\lambda\gtrsim\lambda_{c}^{opt}$, we have combined a perturbation theory with Lagrange multiplier method (LMM) to optimize the prevalence $\rho$. For $\lambda=\lambda_{c}^{opt}+\Delta\lambda$, we have $\mu_{k}=\mu_{k}^{*}+\Delta\mu_{k}$ and $\Theta=\Theta^{*}+\Delta\Theta$, where $\mu_{k}^{*}=\mu k/\left\langle k\right\rangle$, and $\Theta^{*}=1-\frac{\mu}{{\beta\left\langle k\right\rangle}}$ is the solution of self-consistent equation of $\Theta$, Eq.(12) in the main text, under $\mu_{k}=\mu_{k}^{*}$. Since $\Theta>0$ in the region of epidemic spreading, Eq.(12) in the main text can be rewritten as

[TABLE]

Expanding the above equation around ${\left({\mu_{k}^{*},{\Theta^{*}}}\right)}$ to the second-order, it yields

[TABLE]

where $f\buildrel\Delta\over{=}\frac{\beta}{{\left\langle k\right\rangle}}\sum\limits_{k}{\frac{{{k^{2}}P(k)}}{{{\mu_{k}}+\beta k\Theta}}-}1$, and

[TABLE]

Substituting Eq.(S3) into Eq.(A), we obtain

[TABLE]

Using the constraint $\sum\nolimits_{k}P(k)\Delta\mu_{k}=0$ and ignoring the second-order small quantity $\Delta\Theta^{2}\ll\Delta\Theta$, Eq.(S4) becomes

[TABLE]

Around ${\left({\mu_{k}^{*},{\Theta^{*}}}\right)}$, the change $\Delta\rho$ in the prevalence $\rho=\sum\nolimits_{k}P(k)\rho_{k}$ can be expanded in the leading order

[TABLE]

where

[TABLE]

Substituting Eq.(S7) into Eq(S6), we obtain

[TABLE]

In the following we use LMM to minimize $\Delta\rho$ under the constraints $\sum\nolimits_{k}P(k)\Delta\mu_{k}=0$ and Eq.(S5). Note that the first constraint is due to the fixed average curing rate, and the second one is the requirement of the HMF dynamics. Utilizing Eq.(S8) and the two constraints, the Lagrange function can be written as

[TABLE]

where $\tau$ and $\kappa$ are the Lagrange multipliers. Taking the derivative of $\mathcal{L}$ with respect to $\Delta\Theta$ and $\Delta\mu_{k}$, we obtain

[TABLE]

and

[TABLE]

Letting $\frac{{\partial L}}{{\partial\Delta\Theta}}=0$ and $\frac{{\partial L}}{{\partial\Delta\mu_{k}}}=0$, we obtain

[TABLE]

and

[TABLE]

respectively. Substituting Eq.(S13) into the constraint $\sum\nolimits_{k}P(k)\Delta\mu_{k}=0$, we obtain

[TABLE]

Combining Eqs.(S12,S13,S14), we obtain

[TABLE]

Substituting Eq.(S5) and Eq.(S15) into Eq.(S8), we obtain

[TABLE]

Appendix B General Infection Region

For $\lambda$ is larger than but not close to $\lambda_{c}^{opt}$, since the nonlinear character of the model, analytical expression of optimal allocation of $\{\mu_{k}\}$ and the corresponding the minimal $\rho$ is in general impossible. However, with the aid of HMF theory and LMM, the high-dimensional optimization problem can be converted to numerically solving low-dimensional nonlinear equations. We first write a Lagrange function as

[TABLE]

where $\tau$ and $\kappa$ are the Lagrange multipliers. Taking the derivative of $\mathcal{L}$ with respect to $\mu_{k}$ and $\Theta$, we obtain

[TABLE]

and

[TABLE]

Taking the derivative of $\mathcal{L}$ with respect to the Lagrange multipliers $\tau$ and $\kappa$, we obtain the constraint equation Eq.(8) and the self-consistent equation Eq.(12) of $\Theta$ in the main text.

Letting $\partial\mathcal{L}/\partial\mu_{k}=0$, we obtain

[TABLE]

where $k_{c}$ is a threshold degree to guarantee $\mu_{k}>0$ for $k<k_{c}$ and it will be determined later. Substituting Eq.(S20) into Eq.(S19) and letting $\partial\mathcal{L}/\partial\Theta=0$, we obtain

[TABLE]

Combining Eq.(8) in the main text and Eq.(S20), we obtain

[TABLE]

Combining Eq.(12) in the main text and Eq.(S20), we obtain

[TABLE]