Optimal Allocation of Resources for Suppressing Epidemic Spreading on Networks
Hanshuang Chen Guofeng Li Haifeng Zhang Zhonghuai Hou






 chenhshf@ahu.edu.cn





hzhlj@ustc.edu.cn

^1School of Physics and Materials Science, Anhui
University, Hefei, 230601, China 
^2School of Mathematical Science, Anhui University, Hefei, 230601, China 

^3Hefei National Laboratory for Physical Sciences at Microscales
& Department of Chemical Physics, University of
 Science and Technology of China, Hefei, 230026, China






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 λ has a
maximal threshold λ_c^opt=1/⟨ k ⟩
where ⟨ k ⟩ is average degree of the
underlying network. For weak infection region
(λ≳λ_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
(λ>λ_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.


05.10.-a, 64.60.aq, 89.75.Hc Optimal Allocation of Resources for Suppressing Epidemic Spreading on Networks
    Zhonghuai Hou^3
    December 30, 2023
==============================================================================































































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
<cit.>. 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
<cit.> or betweenness centrality
<cit.> to local strategies, like acquaintance
immunization <cit.> and (bias) random walk
immunization <cit.> and to some others
in between <cit.>. Further improvements were done by
graph partitioning <cit.> and the
optimization of the susceptible size <cit.>.
Besides the degree heterogeneity, community structure has also a
major impact on disease immunity
<cit.>. Recently, a message-passing
approach was used to find an optimal set of nodes for immunization
<cit.>. The immunization has been mapped onto the
optimal percolation problem <cit.>. Based on the idea of
explosive percolation, an “explosive immunization" method
 has been proposed <cit.>. 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 <cit.> and
dynamical interplay between awareness and epidemic spreading
<cit.>.


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 <cit.>. 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 <cit.> or the steady-state infection
density <cit.>. 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
𝔸 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 β, and an
infected node i recovers with a nonvanishing curing rate μ_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 μ_i is. Assuming that the total
amount of medicine resource is limited, the average curing rate is
thus fixed, i.e.,

    ⟨μ _i⟩  = μ 10pt  and
     10ptμ _i⩾ 0, 5pt∀ i.

Our goal is to find out an optimal allocation of {μ_i} under
the constraint Eq.(<ref>) so as to minimize the prevalence
ρ, that is the fraction of infected nodes.







In the quenched mean-field (QMF) theory, the probability ρ_i(t)
that node i is infected at time t is described by
N-intertwined equations <cit.>,

    dρ _i(t)/dt =  - μ _iρ _i(t) + β[ 1 - ρ _i(t)]∑_j A_ijρ
    _j(t).

In the steady state, dρ_i(t)/dt=0, ρ_i is determined by a
set of nonlinear equations,

    ρ _i = β∑_j A_ijρ _j/μ _i + β∑_j A_ijρ _j.

One can notice that ρ_i=0 is always a solution of
Eq.(<ref>). This trivial solution corresponds to an absorbing
state with no infective nodes. A nonzero solution ρ_i>0 exists
if the effective infection rate λ=β/μ is larger than
the so-called epidemic threshold λ_c. In this case, the
prevalence ρ = ∑_i ρ _i /N is nonzero
corresponding to an endemic state. By linear stability analysis for
Eq.(<ref>) around ρ_i=0, λ_c is determined by which
the largest eigenvalue of the matrix, -𝕌+β𝔸, is zero, where 𝕌= diag(μ_i) is a diagonal
matrix. For the standard SIS model, μ_i≡μ for all i,
one can immediately obtain the well-known result,
λ_c,QMF^sta=1/Λ_max(𝔸) with the largest
eigenvalue of the adjacency matrix Λ_max(𝔸). 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
𝕌 into two diagonal matrices, 𝕌=𝕌̅+ Δ𝕌, where 𝕌̅= diag{μ
k_i/⟨ k ⟩} with k_i being the degree of
node i and Δ𝕌= diag{Δμ_i}. Since
Tr(𝕌)=Tr(𝕌̅)=Nμ, Δ𝕌 must
satisfy the constraint Tr(Δ𝕌)=0. For the real
symmetric matrix, 𝕌-β𝔸, its largest
eigenvalue Λ_max satisfies the following inequality,

    Λ_max≥v^T (-𝕌+β𝔸)
    v,

where v is a column vector satisfying
v∈ℝ^N and ||v||=1. If we set
v=1/√(N)( 1, ⋯ ,1)^T,
Eq.(<ref>) becomes

    Λ_max≥v^T (-𝕌̅+β𝔸)v- v^T Δ𝕌v
    =-μ+β⟨ k ⟩.

Since Λ_max=0 at the epidemic threshold, Eq.(<ref>)
leads to an upper bound of epidemic threshold,
λ_c≤1/⟨ k ⟩. The condition that
the epidemic threshold equals to the upper bound holds when
v is the eigenvector of 𝕌-β𝔸
corresponding to its largest eigenvalue. If we set
𝕌=𝕌̅ and β=μ/⟨ k
⟩, -𝕌+β𝔸=-μ/⟨ k
⟩𝕃, where 𝕃 is the Laplacian
matrix of the underlying network. It is well-known that the smallest
eigenvalue of 𝕃 is zero and the corresponding eigenvector
is v. Therefore, if the curing rate of each node is
directly proportional to its degree, i.e.,

    μ _i = μ_i^*=μk_i/⟨ k
    ⟩,

the epidemic threshold will be maximized,

    λ_c,QMF^opt=1/⟨ k
    ⟩.

In the QMF theory, the epidemic threshold of the optimal SIS model
is no less than that of the standard SIS model,
λ_c,QMF^opt≥λ_c,QMF^sta, as the lower
bound of Λ_max(𝔸) is ⟨ k
⟩ for any types of networks <cit.>.



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.(<ref>) becomes

    ⟨μ _k⟩  = ∑_k P(k)μ
    _k 1pt  = μ 1pt 1pt 1pt
    1pt 1pt 1pt 1pt 1pt 1pt 1pt 1pt and 1pt 1pt 1pt 1pt 1pt 1pt 1pt 1pt
    1pt 1ptμ _k⩾ 0, 1pt 1pt
    1pt 1pt∀ k 1pt,

where μ_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 <cit.> and its metapopulation version
<cit.> with the curing rate, μ_k∼
k^α, but such a power-law form did not guarantee to be the
optimal one. In <cit.>, 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 <cit.>. 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.(<ref>). The dynamical evolution of ρ_k(t), the
probability of nodes of degree k being infected at time t, reads
<cit.>,

    dρ _k( t )/dt =  - μ _kρ
    _k( t ) + β[ 1 - ρ _k( t )]kΘ(t),

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

    Θ(t)  = ∑_k kP(k)/⟨ k
    ⟩ρ _k(t).

In the steady state, dρ_k(t)/dt=0, Eq.(<ref>) becomes

    ρ _k = β kΘ/μ _k +  β kΘ.

Substituting Eq.(<ref>) into Eq.(<ref>), we obtain a
self-consistent equation of Θ,

    Θ  = ∑_k kP(k)/⟨ k
    ⟩β kΘ/μ _k +  β kΘ.









The epidemic threshold is determined by which the derivation of the
r.h.s of Eq.(<ref>) with respect to Θ at Θ=0
equals to one, leading to

    β _c,HMF = ⟨ k ⟩/∑_k k^2P(k)/μ _k.

For a given P(k), maximizing β_c is equivalent to minimizing
the denominator of the r.h.s of Eq.(<ref>). For this purpose, we
employ Lagrange multiplier method (LMM) to maximize the epidemic
threshold, where the Lagrange function is written as,

    ℒ = ∑_k k^2P(k)/μ _k  +
    τ( ∑_k P(k)μ _k - μ),

where τ is called the Lagrange multiplier. Taking the
derivation of ℒ with respect to μ_k,

    ∂ℒ/∂μ _k =  -
    k^2P(k)/μ _k^2 + τ P(k),

and letting ∂ℒ/
 .
 -∂μ _k = 0 combined with Eq.(8), we arrive at a maximal epidemic threshold

    λ _c,HMF^opt = 1/⟨ k ⟩,

and the corresponding allocation of {μ_k},

    μ _k = μ_k^*=μk/⟨ k ⟩

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,
λ_c,HMF^opt≥λ_c,HMF^sta=⟨ k
⟩/⟨ k^2 ⟩.


For λ larger than but close to λ_c^opt,
λ≳λ_c^opt, we shall combine a perturbation
theory with LMM to optimize the prevalence. To the end, we assume
that for λ = λ_c^opt+ Δλ,
μ_k=μ_k^*+Δμ_k and Θ=Θ^*+ΔΘ,
where Θ^*=1 - μ/β⟨ k ⟩ is the solution of Eq.(<ref>) for μ_k=μ_k^*. Expanding
Eq.(<ref>) around ( μ _k^*,Θ ^*) to
the second order, and then using the constraint ∑_k
P(k)Δμ_k=0 and simultaneously ignoring the second-order
small quantity ΔΘ^2, it yields [See Appendix A for
details]

    ΔΘ  = 1/β ^2⟨ k ⟩∑_k P(k)/kΔμ _k^2
    1pt.

Around ( μ _k^*,Θ ^*), the change Δρ in the prevalence ρ=∑_k P(k)ρ _k can be
written as

    Δρ  =  - Θ ^*/β∑_k
    P(k)/kΔμ _k + ( 1 - Θ ^*)ΔΘ.

Again using LMM to minimize Δρ under the constraints
∑_k P(k)Δμ_k=0 and Eq.(<ref>), we obtain a
minimal ρ=ρ^*+Δρ^opt with ρ ^*
=∑_k P(k)β kΘ ^*/μ _k^* + β
kΘ ^* and

    Δρ^opt   =     - 1/4λ⟨ k
    ⟩ ^2( ⟨k^ - 1⟩ -
    ⟨ k ⟩^ - 1)Δλ ^2
       ≃  - 1/4⟨ k ⟩
    ^3( ⟨k^ - 1⟩  - ⟨
    k ⟩^ - 1)Δλ ^2.

Since ⟨k^ - 1⟩  > ⟨ k
⟩ ^ - 1 for any degree inhomogeneous networks in
terms of Jensen's inequality, Δρ^opt<0 and thus ρ
will be reduced. The corresponding optimal allocation
μ_k=μ_k^*+Δμ_k with

    Δμ _k = μ/2⟨ k ⟩λ( ⟨ k ⟩  - k) Δλ≃μ/2( ⟨ k ⟩-k)Δλ.

This implies that as λ is increased from λ_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, ⟨ k
⟩-k, and the distance of the effective infection rate
to its critical value, Δλ.




For λ is larger than but not close to λ_c^opt,
λ>λ_c^opt, since the nonlinear characteristic of the
model, analytical expression of optimal allocation of {μ_k}
and the corresponding the minimal ρ 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, μ_k satisfies the following
equation,

    μ _k = {√(β kΘ/τ + κβk^2/τ⟨ k ⟩)  - β kΘ >0,    k <k_c 
    
      0,   k ≥ k_c 
    .

where τ and κ are the Lagrange multipliers, and k_c is
a threshold degree to guarantee μ_k>0 for k<k_c and it will be
determined later. Θ, τ and κ are determined by
the following three equations,

    √(βτ⟨ k ⟩/Θ)∑_k = k_min^k_max√(ξ)P(k)  -
    βτ⟨ k ⟩∑_k = k_min^k_maxξ P(k)
     - βκτ∑_k = k_min^k_maxkξ P(k)   -
    κ/⟨ k ⟩Θ∑_k
    = k_c^k_maxkP(k) = 0,


    μ  =  √(βΘ/τ⟨ k
    ⟩)∑_k = k_min^k_ck
    ξ^-1/2P(k) -βΘ∑_k =
    k_min^k_ck P(k),


    Θ= √(βτΘ/⟨ k
    ⟩)∑_k = k_min^k_ck^2
    ξ^-1/2P(k) + 1/⟨ k ⟩∑_k = k_c^k_maxk P(k),

where we have used ξ=k/(⟨ k⟩+κ k).

To numerically solve Θ, τ and κ by
Eqs.(<ref>,<ref>,<ref>), 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 Θ,
τ and κ by Eqs.(<ref>,<ref>,<ref>), and then
test the condition μ_k>0 for all k<k_c by Eq.(<ref>);
(iii) if the condition is not satisfied, k_c will be decreased by
k_c← k_c-1 and return to ii) until the condition
Eq.(<ref>) is fulfilled.





















Figure <ref> shows the optimized results of ρ as a
function λ (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 ⟨ k
⟩=4. For comparison, we also show the results of the
standard SIS model (dotted line) and of the SIS model with the
curing rates μ_i=μ_i^* (dashed line). As expected by the
theoretical prediction, the epidemic threshold of the optimal SIS
model λ_c^opt=1/⟨ k ⟩, which is
significantly larger than that of the standard SIS model, but
coincides with the case of μ_i=μ_i^*. While for
λ>λ_c^opt, the prevalence for μ_i=μ_i^* is
always larger than the optimal choice, and even larger than the
standard SIS model in the strong infection region, indicating that
μ_i=μ_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
ℰ, where ℰ=ρ can be obtained by
numerically iterating Eq.(<ref>). The main steps of SA are as
follows. At beginning, we assign to a given set of {μ_i}
satisfying the constraint Eq.(<ref>) (e.g., μ_i=μ for all
i). Then, we randomly choose two distinct nodes, say i and j,
and try to make the changes μ_i ←μ_i +δ and
μ_j ←μ_j -δ with the standard Metropolis
probability min(1, e^-β_SAΔℰ), where
δ is randomly chosen between -μ_i and μ_i+μ_j to
guarantee the curing rate is always not less than zero. β_SA
is the inverse temperature of SA which slowly increases from
10^-2 to 10^4 via an annealing protocol. Δℰ is the change of the energy function ℰ due the
change of μ_i and μ_j, We tested several different annealing
protocols and we adopted one in which the inverse temperature of SA
β_SA is updated by β_SA← 1.01 β_SA
after each N attempts for updating {μ_i}. The SA results are
also shown in Fig.1 (square dots), which agree with the theoretical
prediction.






In Fig.<ref>, we show the optimal allocation of {μ_k} as
a function of node degree k for several distinct λ 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 λ≳λ_c^opt, μ_k
increases linearly as k with the slope depending on the distance
to the epidemic threshold. The results have been well predicted by
Eq.(<ref>) and Eq.(<ref>). For the region away from the
threshold, μ_k will deviate from linear relation with k. For
sufficiently large λ, μ_k for large k can be less than
that for small k, and even μ_k vanishes when k exceeds a
threshold value, as given by Eq.(<ref>). 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 λ (or the
basic reproduction number R_0=⟨ k ⟩λ). If R_0⩽ 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.(<ref>) or Eq.(<ref>). Once
the maximal epidemic threshold is just across (R_0≳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.(<ref>).
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.(<ref>,<ref>,<ref>). 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
<cit.>.
Moreover, the method presented here could be applied to a number of
other optimization problems, for example, controlling opinion
dynamics in social networks <cit.>. 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).


apsrev


40
natexlab#1#1bibnamefont#1#1bibfnamefont#1#1citenamefont#1#1url<#>1urlprefixURL



[Pastor-Satorras
  et al.(2015)Pastor-Satorras, Castellano, Van Mieghem, and
  Vespignani]RevModPhys.87.925
authorR. Pastor-Satorras,
  authorC. Castellano,
  authorP. Van Mieghem,
  and
  authorA. Vespignani,
  journalRev. Mod. Phys. volume87,
  pages925 (year2015).

[Nowzari et al.(2016)Nowzari, Preciado,
  and Pappas]IEEE.36.26
authorC. Nowzari,
  authorV. M. Preciado,
  and authorG. J.
  Pappas, journalIEEE Control Systems
  volume36, pages26 (year2016).

[Pastor-Satorras and
  Vespignani(2002)]PhysRevE.65.036104
authorR. Pastor-Satorras
  and
  authorA. Vespignani,
  journalPhys. Rev. E volume65,
  pages036104 (year2002).

[Holme et al.(2002)Holme, Kim, Yoon, and
  Han]PhysRevE.65.056109
authorP. Holme,
  authorB. J. Kim,
  authorC. N. Yoon, and
  authorS. K. Han,
  journalPhys. Rev. E volume65,
  pages056109 (year2002).

[Cohen et al.(2003)Cohen, Havlin, and ben
  Avraham]PhysRevLett.91.247901
authorR. Cohen,
  authorS. Havlin, and
  authorD. ben Avraham,
  journalPhys. Rev. Lett. volume91,
  pages247901 (year2003).

[Holme(2004)]EPL.68.908
authorP. Holme,
  journalEurophys. Lett. volume68,
  pages908 (year2004).

[Stauffer and Barbosa(2006)]PhysRevE.74.056105
authorA. O. Stauffer and
  authorV. C. Barbosa,
  journalPhys. Rev. E volume74,
  pages056105 (year2006).

[Gomez-Gardenes
  et al.(2006)Gomez-Gardenes, Echenique, and
  Moreno]EPJB49.259
authorJ. Gomez-Gardenes,
  authorP. Echenique,
  and authorY. Moreno,
  journalEur. Phys. J. B volume49,
  pages259 (year2006).

[Chen et al.(2008)Chen, Paul, Havlin,
  Liljeros, and Stanley]PhysRevLett.101.058701
authorY. Chen,
  authorG. Paul,
  authorS. Havlin,
  authorF. Liljeros, and
  authorH. E. Stanley,
  journalPhys. Rev. Lett. volume101,
  pages058701 (year2008).

[Schneider et al.(2011)Schneider,
  Mihaljev, Havlin, and Herrmann]PhysRevE.84.061911
authorC. M. Schneider,
  authorT. Mihaljev,
  authorS. Havlin, and
  authorH. J. Herrmann,
  journalPhys. Rev. E volume84,
  pages061911 (year2011).

[Masuda(2009)]NJP11.123018
authorN. Masuda, journalNew
  J. Phys. volume11, pages123018
  (year2009).

[Salathé and Jones(2010)]PlosCP6.e1000736
authorM. Salathé and
  authorJ. H. Jones,
  journalPLoS Comput. Biol. volume6,
  pagese1000736 (year2010).

[Altarelli et al.(2014)Altarelli,
  Braunstein, Dall'Asta, Wakeling, and Zecchina]PhysRevX.4.021024
authorF. Altarelli,
  authorA. Braunstein,
  authorL. Dall'Asta,
  authorJ. R. Wakeling,
  and authorR. Zecchina,
  journalPhys. Rev. X volume4,
  pages021024 (year2014).

[Morone and Makse(2015)]Morone2015
authorF. Morone and
  authorH. A. Makse,
  journalNature (London) volume524,
  pages65 (year2015).

[Clusella et al.(2016)Clusella,
  Grassberger, Pérez-Reche, and Politi]PhysRevLett.117.208301
authorP. Clusella,
  authorP. Grassberger,
  authorF. J. Pérez-Reche,
  and authorA. Politi,
  journalPhys. Rev. Lett. volume117,
  pages208301 (year2016).

[Gross et al.(2006)Gross, D'Lima, and
  Blasius]PhysRevLett.96.208701
authorT. Gross,
  authorC. J. D. D'Lima,
  and authorB. Blasius,
  journalPhys. Rev. Lett. volume96,
  pages208701 (year2006).

[Granell et al.(2013)Granell, Gómez,
  and Arenas]PhysRevLett.111.128701
authorC. Granell,
  authorS. Gómez, and
  authorA. Arenas,
  journalPhys. Rev. Lett. volume111,
  pages128701 (year2013).

[Prakash et al.(2013)Prakash, Adamic,
  Iwashnya, Tong, and Faloutsos]Prakash2013
authorB. A. Prakash,
  authorL. Adamic,
  authorT. Iwashnya,
  authorH. Tong, and
  authorC. Faloutsos,
  journalProc. SIAM Int. Conf. Data Mining, Austin, TX p.
  pages659C667 (year2013).

[Wan et al.(2008)Wan, Roy, and
  Saberi]Wan2008
authorY. Wan,
  authorS. Roy, and
  authorA. Saberi,
  journalSyst. Biol. IET volume2,
  pages184 (year2008).

[Preciado et al.(2013)Preciado, Zargham,
  Enyioha, Jadbabaie, and Pappas]Preciado2013
authorV. M. Preciado,
  authorM. Zargham,
  authorC. Enyioha,
  authorA. Jadbabaie,
  and authorG. J.
  Pappas, journalProc. IEEE Conf. Decision Control,
  Florence, Italy pp. pages7486–7491 (year2013).

[Gourdin et al.(2011)Gourdin, Omic, and
  Mieghem]Gourdin2011
authorE. Gourdin,
  authorJ. Omic, and
  authorP. V. Mieghem,
  journalProc. 8th Int. Workshop Design Reliable Communication
  Networks pp. pages86–93 (year2011).

[Y. Wang()]Wang2003
authorY. Wang et al.,
  journal22nd International Symposium on Reliable Distributed
  Systems (SRDS03) (IEEE Computer Society, 2003), p. 25.

[Mieghem et al.(2009)Mieghem, Omic, and
  Kooij]Mieghem2009
authorP. V. Mieghem,
  authorJ. Omic, and
  authorR. Kooij,
  journalIEEE ACM Trans. Netw. volume17,
  pages1 (year2009).

[Gómez et al.(2010)Gómez, Arenas,
  Borge-Holthoefer, Meloni, and Moreno]Gomez2010
authorS. Gómez,
  authorA. Arenas,
  authorJ. Borge-Holthoefer,
  authorS. Meloni, and
  authorY. Moreno,
  journalEurophys. Lett. volume89,
  pages38009 (year2010).

[Mieghem(2011)]Mieghem_Book2011
authorP. V. Mieghem,
  titleGraph Spectra for Complex Networks
  (publisherCambridge University Press, year2011).

[Dezs ő and
  Barabási(2002)]PhysRevE.65.055103
authorZ. Dezs ő and authorA.-L.
  Barabási, journalPhys. Rev. E
  volume65, pages055103
  (year2002).

[Shen et al.(2012)Shen, Chen, and
  Hou]PhysRevE.86.036114
authorC. Shen,
  authorH. Chen, and
  authorZ. Hou,
  journalPhys. Rev. E volume86,
  pages036114 (year2012).

[Borgs et al.(2010)Borgs, Chayes, Ganesh,
  and Saberi]Borgs2010
authorC. Borgs,
  authorJ. Chayes,
  authorA. Ganesh, and
  authorA. Saberi,
  journalRandom Struct. Alg. volume37,
  pages204 (year2010).

[Chung et al.(2009)Chung, Horn, and
  Tsiatas]Chung2009
authorF. Chung,
  authorP. Horn, and
  authorA. Tsiatas,
  journalInternet Math. volume6,
  pages237 (year2009).

[Pastor-Satorras and Vespignani(2001)]PRL01003200
authorR. Pastor-Satorras
  and
  authorA. Vespignani,
  journalPhys. Rev. Lett. volume86,
  pages3200 (year2001).


[Eames and Keeling(2002)]Eames2002
authorK. T. D. Eames and
  authorM. J. Keeling,
  journalProc. Natl. Acad. Sci. U.S.A.
  volume99, pages13330 (year2002).

[Gleeson(2011)]PhysRevLett.107.068701
authorJ. P. Gleeson,
  journalPhys. Rev. Lett. volume107,
  pages068701 (year2011).

[Boguñá et al.(2013)Boguñá,
  Castellano, and Pastor-Satorras]PhysRevLett.111.068701
authorM. Boguñá,
  authorC. Castellano,
  and
  authorR. Pastor-Satorras,
  journalPhys. Rev. Lett. volume111,
  pages068701 (year2013).

[Mata et al.(2014)Mata, Ferreira, and
  Ferreira]Mata2014
authorA. S. Mata,
  authorR. S. Ferreira,
  and authorS. C.
  Ferreira, journalNew J. Phys.
  volume16, pages053006
  (year2014).

[Kiss et al.(2015)Kiss, Röst, and
  Vizi]PhysRevLett.115.078701
authorI. Z. Kiss,
  authorG. Röst, and
  authorZ. Vizi,
  journalPhys. Rev. Lett. volume115,
  pages078701 (year2015).

[Cai et al.(2016)Cai, Wu, Chen, Holme,
  and Guan]PhysRevLett.116.258301
authorC.-R. Cai,
  authorZ.-X. Wu,
  authorM. Z. Q. Chen,
  authorP. Holme, and
  authorJ.-Y. Guan,
  journalPhys. Rev. Lett. volume116,
  pages258301 (year2016).

[Cator and Van Mieghem(2012)]PhysRevE.85.056111
authorE. Cator and
  authorP. Van Mieghem,
  journalPhys. Rev. E volume85,
  pages056111 (year2012).

[Mata and Ferreira(2013)]Mata2103
authorA. S. Mata and
  authorS. C. Ferreira,
  journalEurophys. Lett. volume103,
  pages48003 (year2013).

[I. Z. Kiss(2015)]Kiss2015
author I. Z. Kiss et al.,
  journalJ. Math. Biol. volume70,
  pages437 (year2015).

[Castellano et al.(2009)Castellano,
  Fortunato, and Loreto]RMP09000591
authorC. Castellano,
  authorS. Fortunato,
  and authorV. Loreto,
  journalRev. Mod. Phys. volume81,
  pages591 (year2009).






§ WEAK INFECTION REGION

For λ larger than but close to λ_c^opt, λ≳λ_c^opt, we have combined a perturbation theory
with Lagrange multiplier method (LMM) to optimize the prevalence
ρ. For λ = λ_c^opt+ Δλ, we have
μ_k=μ_k^*+Δμ_k and Θ=Θ^*+ΔΘ,
where μ_k^*=μ k/⟨ k ⟩, and Θ^*=1
- μ/β⟨ k ⟩ is the
solution of self-consistent equation of Θ, Eq.(<ref>) in
the main text, under μ_k=μ_k^*. Since Θ>0 in the region
of epidemic spreading, Eq.(<ref>) in the main text can be
rewritten as

    S1β/⟨ k ⟩∑_k
    k^2P(k)/μ _k + β kΘ=1.

Expanding the above equation around ( μ _k^*,Θ ^*) to the second-order, it yields

    S2   ∑_k . ∂ f/∂μ _k|_( μ _k^*,Θ ^*)Δμ _k +
    . ∂ f/∂Θ|_(
    μ _k^*,Θ ^*)ΔΘ
       +
    1/2∑_k . ∑_k'∂ ^2f/∂μ _k∂μ _k'|_( μ _k^*,Θ ^*)Δμ
    _kΔμ _k' 1pt 1pt 1pt
        + ∑_k . ∂
    ^2f/∂μ _k∂Θ|_( μ
    _k^*,Θ ^*)Δμ _kΔΘ  + .
    1/2∂ ^2f/∂Θ∂Θ|_( μ _k^*,Θ ^*)ΔΘ^2
    =0,

where f Δ = β/⟨ k
⟩∑_k k^2P(k)/μ _k + β
kΘ -  1, and

    S3   . ∂ f/∂μ _k|_( μ _k^*,Θ ^*) =   - P(k)/β⟨ k ⟩
       . ∂ f/∂Θ|_( μ _k^*,Θ ^*) =   - 1
       . ∂ ^2f/∂μ _k∂μ _k'|_( μ _k^*,Θ ^*) = δ _kk'2P(k)/β ^2⟨ k ⟩ k
       . ∂ ^2f/∂μ _k∂Θ|_( μ _k^*,Θ ^*) = 2P(k)/β⟨ k ⟩
       . ∂ ^2f/∂Θ∂Θ|_( μ _k^*,Θ ^*) = 2.

Substituting Eq.(<ref>) into Eq.(<ref>), we obtain

    S4    - 1/β⟨ k ⟩∑_k
    P(k)Δμ _k - ΔΘ  + 1/β
    ^2⟨ k ⟩∑_k P(k)/kΔμ _k^2 1pt
       + 2/β⟨ k ⟩∑_k
    P(k)Δμ _kΔΘ  + ΔΘ ^2 = 0.

Using the constraint ∑_k P(k)Δμ_k=0 and
ignoring the second-order small quantity ΔΘ^2≪ΔΘ, Eq.(<ref>) becomes

    S5ΔΘ  = 1/β ^2⟨ k ⟩∑_k P(k)/kΔμ _k^2 1pt.

Around ( μ _k^*,Θ ^*), the change Δρ in the prevalence ρ=∑_k P(k)ρ _k can be
expanded in the leading order

    S6Δρ  = ∑_k . ∂ρ/∂μ _k|_( μ _k^*,Θ ^*)Δμ _k + . ∂ρ/∂Θ|_( μ _k^*,Θ ^*)ΔΘ,

where

    S7   . ∂ρ/∂μ _k|_( μ _k^*,Θ ^*) =  - P(k)Θ ^*/β k, 
       . ∂ρ/∂Θ|_( μ _k^*,Θ ^*) = 1 - Θ ^*.

Substituting Eq.(<ref>) into Eq(<ref>), we obtain

    S8Δρ  =  - Θ ^*/β∑_k
    P(k)/kΔμ _k + ( 1 - Θ ^*)ΔΘ.

In the following we use LMM to minimize Δρ under the
constraints ∑_k P(k)Δμ_k=0 and Eq.(<ref>).
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.(<ref>) and the two constraints, the Lagrange
function can be written as

    S9ℒ =     - Θ ^*/β∑_k
    P(k)/kΔμ _k + ( 1 - Θ ^*)ΔΘ
       + τ(  - ΔΘ  +
    1/β ^2⟨ k ⟩∑_k
    P(k)/kΔμ _k^2) 1pt  + κ∑_k P(k)Δμ _k,

where τ and κ are the Lagrange multipliers. Taking the
derivative of ℒ with respect to ΔΘ and
Δμ_k, we obtain

    S10∂ℒ/∂ΔΘ = (
    1 - Θ ^*) - τ,

and

    S11∂ℒ/∂Δμ _k =  -
    Θ ^*/βP(k)/k + τ1/β
    ^2⟨ k ⟩2P(k)/kΔμ _k +
    κ P(k).

Letting ∂ L/∂ΔΘ =0 and
∂ L/∂Δμ_k  =0, we obtain

    S12τ  = 1 - Θ ^*,

and

    S13
     - Θ ^*/β k + 2τ/β
    ^2k⟨ k ⟩Δμ _k + κ  = 0,

respectively. Substituting Eq.(<ref>) into the constraint
∑_k P(k)Δμ_k=0, we obtain

    S14κ  = Θ ^*/β⟨ k ⟩.

Combining Eqs.(<ref>,<ref>,<ref>), we obtain

    S15Δμ _k = μ/2⟨ k ⟩λ( ⟨ k ⟩  - k) Δλ≃μ/2( ⟨ k ⟩  - k)Δλ.

Substituting Eq.(<ref>) and Eq.(<ref>) into Eq.(<ref>),
we obtain

    S16Δρ^opt  = - 1/4λ⟨ k
    ⟩ ^2( ⟨k^ - 1⟩ -
    ⟨ k ⟩^ - 1)Δλ ^2
    ≃  - 1/4⟨ k ⟩
    ^3( ⟨k^ - 1⟩  - ⟨
    k ⟩^ - 1)Δλ ^2.





§ GENERAL INFECTION REGION

For λ is larger than but not close to λ_c^opt,
since the nonlinear character of the model, analytical expression of
optimal allocation of {μ_k} and the corresponding the minimal
ρ 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

    S17ℒ=     ∑_k P(k)β kΘ/μ
    _k + β kΘ + τ( ∑_k P(k)μ _k
    - μ) 
       + κ( ∑_k
    kP(k)/⟨ k ⟩β k
    Θ/μ _k + β kΘ - Θ),

where τ and κ are the Lagrange multipliers. Taking the
derivative of ℒ with respect to μ_k and Θ,
we obtain

    S18∂ℒ/∂μ _k =  -
    P(k)β kΘ/( μ _k + β kΘ)^2 + τ P(k)- κkP(k)/⟨ k
    ⟩β k Θ/( μ _k + β
    k )^2,

and

    S19∂ℒ/∂Θ =     ∑_k
    β kP(k)/μ _k + β kΘ - ∑_k
    β ^2k^2P(k)Θ/( μ _k + β
    kΘ)^2
        - κ∑_k
    kP(k)/⟨ k ⟩β
    ^2k^2Θ/( μ _k + β kΘ)^2.

Taking the derivative of ℒ with respect to the Lagrange
multipliers τ and κ, we obtain the constraint equation
Eq.(<ref>) and the self-consistent equation Eq.(<ref>) of
Θ in the main text.

Letting ∂ℒ/∂μ_k=0, we obtain

    S20μ _k = {√(β kΘ/τ + κβk^2/τ⟨ k ⟩)  - β kΘ >0,    k <k_c 
    
      0,   k ≥ k_c 
    .

where k_c is a threshold degree to guarantee μ_k>0 for k<k_c
and it will be determined later. Substituting Eq.(<ref>) into
Eq.(<ref>) and letting ∂ℒ/∂Θ=0, we obtain

    S21√(βτ⟨ k ⟩/Θ)∑_k = k_min^k_cP(k)√(k)/√(⟨ k ⟩  + κ k)  - βτ⟨ k ⟩∑_k = k_min^k_cP(k)k/⟨ k ⟩  + κ k
    - βκτ∑_k = k_min^k_cP(k)k^2/⟨ k ⟩
     + κ k  - κ/⟨ k
    ⟩Θ∑_k = k_c^k_maxkP(k)  = 0.

Combining Eq.(<ref>) in the main text and Eq.(<ref>), we
obtain

    S22μ  =     √(βΘ/τ⟨ k
    ⟩)∑_k = k_min^k_cP(k)√(k)( √(⟨ k ⟩  + κ k)) 
       -βΘ∑_k = k_min^k_cP(k)k.

Combining Eq.(<ref>) in the main text and Eq.(<ref>), we
obtain

    S23Θ= √(βτΘ/⟨ k
    ⟩)∑_k = k_min^k_cP(k)k^3/2/√(⟨ k ⟩  +
    κ k)  + 1/⟨ k ⟩∑_k = k_c^k_maxkP(k)  .




