Containing epidemic outbreaks by message-passing techniques

The problem of targeted network immunization can be deﬁned as the one of ﬁnding a subset of nodes in a network to immunize or vaccinate in order to minimize a tradeoﬀ between the cost of vaccination and the ﬁnal (stationary) expected infection under a given epidemic model. Although computing the expected infection is a hard computational problem, simple and eﬃcient mean-ﬁeld approximations have been put forward in the literature in recent years. The optimization problem can be recast into a constrained one in which the constraints enforce local mean-ﬁeld equations describing the average stationary state of the epidemic process. For a wide class of epidemic models, including the susceptible-infected-removed and the susceptible-infected-susceptible models, we deﬁne a message-passing approach to network immunization that allows us to study the statistical properties of epidemic outbreaks in the presence of immunized nodes as well as to ﬁnd (nearly) optimal immunization sets for a given choice of parameters and costs. The algorithm scales linearly with the size of the graph and it can be made eﬃcient even on large networks. We compare its performance with topologically based heuristics, greedy methods, and simulated annealing.


I. INTRODUCTION
One of the key questions of computational epidemiology is how best to distribute limited resources of treatment and vaccination so that they will be most effective in suppressing or reducing outbreaks of disease.This problem is heightened by the entangled networks of interactions via which diseases can spread: in a large complex network, contact with a high-degree hub can see a virus spread rapidly throughout the population even if the probability of transmission from an individual contact is low [1,2].Early works on network immunization drew attention to the differences between random immunization and targeted immunization strategies [3][4][5][6].A simple random immunization strategy can consist in fixing a fraction or a density of immunized nodes and averaging the outcome of the epidemic process over all possible realizations of the immunization set.On the contrary, targeted immunization strategies correlate the choice of immunized nodes with some topological feature, such as the degree or other centrality measures.This can be experimentally shown to have some positive effect in reducing the spread of diseases [3][4][5][6].Most topologically-based algorithms for immunization follow an incremental procedure, in which the set of immunized nodes is initially empty then it is progressively filled adding one by one the nodes that are most relevant with respect to a particular topological metric.Despite the computational cost, recalculating the topological metric after each immunization step (i.e. after removing the immunized node from the graph) usually provides much better results than computing it only once on the original graph [4].Further improvements were obtained by means of more complex immunization strategies, based on graph partitioning [7,8] and on the optimization of the susceptible size [9].Beside the heterogeneity of contacts, also clustering, community structure and modularity have a major impact on disease dynamics [10,11], therefore the same immunization strategy can produce contrasting results on networks with different topological features.This is a consequence of the fact that topological heuristic methods neglect important features of the spreading rule and most common metrics used to measure their effectiveness, such as the largest connected component or the largest non-immune cluster size, are proxies that may not reflect the true susceptibility to an epidemic.These techniques also neglect the cost of vaccination, which may vary widely depending on the chosen target.
To overcome these limitations, several authors tried to quantify more explicitly the effects of immunization strategies on the outbreak dynamics [12][13][14][15][16][17] and network immunization was mathematically formulated as a proper optimization problem, that can be proven to be NP-hard in a plethora of different variants [16,[18][19][20][21][22][23][24][25].Standard optimization techniques such as Monte-Carlo (MC) methods or integer/linear programming are computationally very expensive and may take a prohibitive amount of time to reach reasonably good results even on relatively small networks.On the other hand, the greedy optimization strategies usually proposed are guaranteed to approximate the optimal result by a constant factor only in some fortunate case [16,18,25].
Recent progress in combinatorial optimization have shown that algorithms based on the message-passing principle, and developed using methods from the statistical physics of disordered systems, outperform in many cases both greedy algorithms and simulated annealing, even in complex optimization problem involving stochastic parameters [26,27] and dynamical rules [28,29].In many cases in which MC algorithms get trapped in local minima of the (free-)energy function, message-passing algorithms can find considerably better results.The remarkable performances of these algorithms are combined with considerably good computation time scaling properties.While on a wide variety of optimization problems the computational complexity of simulated annealing scales exponentially with the system size, message-passing algorithms typically require a time that scales roughly linearly with the number of messages (i.e. the number of edges).
In this paper we show that, under some approximations, network immunization can be written as a constrained optimization problem, in which the constraints are fixed-point equations for some local (node or edge) variables describing the average stationary state of the dynamics.These constraints and a suitably defined objective (energy) function are then used to derive a message-passing approach to the optimization problem, and to design efficient algorithms on large networks.We apply this method to find optimal immunization strategies for both susceptibleinfected-recovered (SIR) and susceptible-infected-susceptible (SIS) models.
In Section II we recall the main ideas and formulas of mean-field methods in epidemic models, that are usually used to estimate the average stationary properties of an epidemic outbreak.The optimal immunization problem is introduced in Section III, opportunely defined in terms of mean-field quantities.Section IV is devoted to the definition of the message-passing approach and the derivation of the corresponding Belief-Propagation (BP) and Max-Sum (MS) equations.In Section V, we use BP to understand in detail the immunization properties on the prototypical case of random regular graphs.The comparison with other optimization methods on more general graphs is discussed in Section VI.

II. MEAN-FIELD METHODS IN EPIDEMIC MODELS
Over the years, a large number of stochastic epidemic models have been introduced, with the aim of addressing some specific features of different diseases [30,31].In the most simple model, the epidemic spreading induces in the nodes irreversible stochastic transitions from a susceptible state to an infected one.Infected individuals can recover either returning to the susceptible state or becoming permanently resistant to the disease.One can then increase the complexity of the stochastic model introducing other intermediate states, or compartments, such as exposure and latency.In the following, we discuss the most basic models of epidemic spreading, providing for each of them a set of approximated equations of mean-field type valid on very general graph structures.Their solution describes the statistical properties of the stationary state corresponding to a given set of initial conditions and external parameters.In addition, such equations allow to measure the level of infection once a configuration of initially immunized nodes is chosen.

A. Irreversible Epidemic Processes
The susceptible-infected-recovered (SIR) model was formulated by Kermack and McKendrick [32] to describe the irreversible propagation through a population of individuals of an infectious disease, such as measles, mumps, or cholera.The SIR stochastic dynamics is defined over a graph G = (V, E), representing the contact network of a set V of individuals.At any given time step t (e.g. a day), a node i can be in one of three states: susceptible (S), infected (I), and recovered/removed (R).The state of node i at time t is represented by a variable x t i ∈ {S, I, R}.We assume that each node i is initially infected with probability q i ∈ [0, 1] (independent of the other nodes).At each time step, an infected node i can first spread the disease to each susceptible neighbor j with given probability T ij ∈ (0, 1], then recover with probability r i .Once recovered, individuals do not get sick anymore (they are effectively removed from the graph).The probability that an infected node i directly transmits the disease to j before i recovers is given by It is thus possible to construct a completely static representation of the process that maps the final state onto the outcome of a bond percolation process [33].This relationship can be made mathematically clear as follows.Let us consider a tree-like graph and define m ij to be the probability that node i is eventually infected when considering the graph obtained in the absence of the neighboring node j.Exploiting the factorization of probabilities on the sub-branches of the tree emerging from i, the quantity m ij satisfies the equation where ∂i denotes the set of neighbors of i.Since infected nodes eventually recover, in the final state nodes can only be either in the susceptible state or in the recovered one.From the knowledge of the conditional marginals m ij , one can compute the probability m i that a node i is eventually infected, i.e. the probability that i is recovered in the final state, Although ( 2)-( 3) are exact only on trees, they have been successfully applied to study the SIR model also on general random graphs [33][34][35].A comparison between the solutions of these equations and the results of simulations of the SIR stochastic process is shown in Figure 1 for a random regular graph (RRG) of N = 10 3 nodes and degree K = 4.For simplicity we considered uniform self-infection probabilities q i = q, ∀i ∈ V and uniform transmission probabilities p ij = p, ∀(i, j) ∈ E. In the SIR stochastic process, we defined a measure of the "outbreak" size as the average fraction f of nodes that have been infected during the epidemic spreading.Since all infected nodes eventually recover, this metric can be also defined as where Pr [x ∞ i = R] is the probability that node i is eventually recovered.In Fig. 1, f is plotted as function of the transmission probability p for q = 0.1, 0.01, 0.001.The results, obtained averaging over 10 4 realization of the stochastic process, are reported as black circles, while the symmetric bars indicate the fluctuations around the average behavior.The average behavior can also be computed from the solutions of ( 2)-( 3), exploiting the fact that, in the tree-like approximation Pr [x ∞ i = R] m i .The results are reported as a red full line in Fig. 1.The agreement between the mean-field theory represented by ( 2)-(3) and the simulations is very good for sufficiently large values of q, then it deteriorates for q of the order of 1/N and large values of p.The reason for such a discrepancy is that Eqs.(2) are correct on tree-like structures, i.e. when the disease transmission events to one node coming from two neighbors are not correlated.The "decorrelation" assumption is not correct when the actual number of sources of spontaneous infection is very small.This is obvious in the case of a unique source of infection: the contagion path has the same source, hence the infection of a node due to disease transmission from her neighbors is a highly correlated process, that is not well captured by (2)-(3).More precisely the solution to Eqs. ( 2)-(3) gives an upper bound for the real probability to be infected [35].In the limit of infinitely large networks, this approach is expected to provide a correct description of the average final state of the system for any finite value of q and p.
Equations ( 2)-( 3) can be easily modified to include immunization of nodes.By considering a set of binary variable s i ∈ {0, 1}, in which s i = 1 if node i is immune to the disease, we get ∀i ∈ V and Given a configuration s = {s 1 , . . ., s N } of immune nodes, that we call immunization set, the solution of (5a)-(5b) provides a measure of the corresponding epidemic outbreak.It is possible to show that for a given set of parameters {q i } and {p ij }, the solution of (2)-( 3) is unique, therefore each configuration of immune nodes s corresponds to a unique solution of the equations (5a)-(5b).This property will be crucial for the validity of the optimization method developed in this work.

B. Reversible Epidemic Processes
The susceptible-infected-susceptible (SIS) model is the prototype of reversible models of epidemic spreading, in which after recovery a node is again susceptible of being infected [30,31].The state of node i at time t is now represented by a binary variable x t i ∈ {S, I}.At each time step, an infected node i can transmit the disease to each of its susceptible neighbors j with probability p ij , while it recovers with rate r i (becoming susceptible again).The stochastic process admits an absorbing state in which all nodes are susceptible and the disease has disappeared from the population.When the transmission probabilities are sufficiently large, an active stationary state also exists, that is metastable and attractive for the dynamical process.Although in any finite population a fluctuation will eventually bring the system into the absorbing state, the lifetime of the metastable endemic state scales with the size of the graph in such a way that an absorbing phase transition as function of the transmission probabilities is expected to occur in the thermodynamic limit [36][37][38][39][40].The critical threshold usually depends on the parameters of the dynamical process as well as on the topological structure of the underlying interaction graph.A variant of the SIS model with spontaneous self-infection was recently introduced [41,42] in order to simplify the numerical and mathematical analysis of the model.The presence of spontaneous self-infection destroys the absorbing state, whereas the metastable state becomes the (unique) stationary state of the dynamics.On the other hand, for a given small self-infection probability, the dynamics shows a clear boundary between low infection region and a region of global spreading as function of the transmission probabilities.By scaling down self-infection, one can extrapolate information on the epidemic phase transition occurring for zero self-infection, avoiding the problems associated with the existence of an absorbing state.
The SIS model on a given graph with N nodes is a Markov chain with 2 N states, whose stationary probability distribution cannot be explicitly computed for large systems.A simple mean-field approximation that turned out to provide a good qualitative and quantitive description of the stationary state of the SIS model is obtained replacing the exact probability distribution by a product measure over the nodes of the graph [38][39][40][43][44][45].This factorization, also known as the N -intertwined model, leads to a set of "quenched" mean-field equations for the evolution of single-node infection probabilities m i with i ∈ V .Notice the difference between the SIR and SIS cases: while in the SIR model m i indicates the (mean-field) probability that node i is eventually infected before the final state, in the SIS model it represents the (mean-field) probability that i is infected in the stationary state of the dynamics.In the "quenched" mean-field approximation, the infection probability of node i at time t satisfies the equation where p ji is the transmission probability from j to i, q i is the spontaneous self-infection probability.In the stationary state, the mean-field variables {m i } i∈V are given by the solution of the fixed-point equations where s i ∈ {0, 1} says whether node i is immune or not.As for the SIS model, it is possible to show that the mean-field quantity m i gives an upper bound for the real value of the infection probability of node i in the stationary state [38][39][40].The approximation can be improved considering second-order quantities, i.e. deriving closed equations for single-point marginals and pair-correlations, but the actual form of these equations is not unique and depends on the moment closure approximation adopted [38][39][40].Nevertheless, in most cases equations (7) already provide a very good description of the stationary state of the SIS stochastic process.
A measure of the outbreak size of the epidemics is given by the average fraction f of infected nodes in the active stationary state.If the stationary state is infinitely long-lived, f can be operatively defined as where 1[•] is equal to 1 when the argument is verified and 0 otherwise.In practice, simulations are performed for a finite time, that is chosen to be much longer than the time necessary to converge to the stationary state.Figure 2 displays the behavior of f as function of the transmission probabilities p, for different values of the self-infection probability q, on a random regular graph of N = 10 3 nodes and degree K = 4.The dashed line is the same quantity computed as f i∈V m i , i.e. approximating the probability that node i is infected in the stationary state by the mean-field one m i obtained solving (7).The agreement between mean-field predictions and results of the simulations is very good in almost all regimes of the parameters.

C. Continuous-Time Processes
Both SIR and SIS epidemic models are often formulated in continuous time, where spontaneous self-infection probabilities and transmission probabilities are replaced by rates.While the time-dependent evolution of the dynamical processes can be quite different from that of their discrete-time counterpart, the stationary behavior at long-time can be similarly described in terms of percolation-like equations.In continuous time, transmission events can be modeled as independent Poisson processes, therefore the probability that, in a sufficiently small interval of time, a node is infected from the neighbors is proportional to the sum of the independent probabilities of the individual transmission events.Once {q i } i∈V and {p ij } (i,j)∈E are intended as probability rates, this reduces to replace the term 1 − j (1 − p ji m j ) with j p ji m j in both (5a)-(5b) and (7).

III. THE OPTIMAL IMMUNIZATION PROBLEM
Preventing or eradicating diseases entails a trade-off between the costs of treating and hospitalizing infected individuals and the cost of distributing vaccines/drugs.When the contact network is known and these costs can be estimated, one can devise an optimization problem, in which the optimal immunization set is the configuration of immunized nodes that minimizes a properly defined energy function.In the SIR stochastic process we consider the following energy function in which c i ∈ R is the cost of immunization of node i, i is a loss, i.e. the cost associated with the infection of node i, while Pr [x ∞ i = R|s] is the probability that node i eventually recovers (i.e. the probability that the node has been infected during the epidemic spreading) given the configuration s of immunized nodes.Estimating the probability Pr [x ∞ i = R|s] from the simulations of the stochastic SIR process is very cumbersome, making the optimization problem practically unsolvable for sufficiently large graphs.In fact, it is known that finding the probability that a node becomes infected in the course of an epidemic in SIR and related models is an NP-hard problem [46].The energy function for the optimal immunization problem for the SIS model can be defined in a similar way to be Likewise the SIR case, computing E SIS (s) from direct simulations of the stochastic dynamics is a very demanding task.
The problem of finding the configuration of seeds that minimizes an energy function associated with a given dynamical model has been investigated by several authors in the recent past.Even for simple stochastic propagation models, such as the independent cascade model [19], several versions of this optimization problem have been introduced and shown to be NP-hard.To the best of our knowledge, the methods proposed to solve this class of problems are mostly based on greedy approaches, that scale reasonably well with the system size but give solutions that only in very peculiar cases provide good approximations of the real optima.
In this work we employ a simplified approach, based on the mean-field representation of the dynamical process.We replace the energy functions E SIR and E SIS with a unique energy function defined as in which ∀i ∈ V , m i is the solution of (5a)-( 5b) and ( 7) for the SIR and SIS models respectively.This energy can now be computed easily for any configuration s of immunized nodes, solving the corresponding set of self-consistent equations for {m i }.Once we know how to compute the energy function, the optimization problem can be studied with different methods.The simplest optimization method is a greedy algorithm that resembles the incremental procedure used in topologically-based heuristics.Starting from an empty set of immunized nodes, the greedy principle imposes to select the node whose immunization causes the largest drop in the energy function.Once the node is included in the set of immunized nodes, a second node is chosen following the same greedy principle.One-by-one all nodes are added to the immunization set.Under this criterion, the best configuration of immunized nodes is the one that realizes the lowest energy.The greedy algorithm is fast but it presents several drawbacks: first, there is no guarantee that the best solution found is the optimal one; moreover, we cannot stop the greedy procedure before almost all nodes are added to the immunization set, because the energy function can display several local minima during this process.Another famous optimization method is simulated annealing, in which a MC algorithm is used to find the configuration s minimizing the energy E (s).For a sufficiently slow annealing schedule, the algorithm is guaranteed of finding the minimum, i.e. the optimal immunization set.However, this can be computationally unfeasible in large graphs, since MC algorithms get trapped in local minima of the (free-)energy function and the convergence to the global one requires a time that can scale exponentially with the system size.In the next section, we will develop a message-passing approach to the optimal immunization problem.Message-passing algorithms can find the global optimum or, in general, considerably better results than simulated annealing in a running time that scales only linearly with the number of messages (i.e. the number of edges).The implementation details of the greedy algorithm and simulated annealing are described in Appendix A.

IV. BELIEF-PROPAGATION APPROACH TO THE OPTIMAL IMMUNIZATION PROBLEM
Since each configuration of immunized nodes corresponds to a unique value of the energy function, we can define a probability function Q(s) by associating to each configuration s a Boltzmann weight e −βE (s,m) and tracing over the variables m with the constraint that they are a solution of the corresponding percolation-like equations.It follows that the equations in (5a) and in (7) become a set of hard constraints for the auxiliary variables m that must be included in the optimization problem in addition to the energetic terms.As for standard constraint-satisfaction problems over Factor graph representation of the BP equations for the optimal immunization problem in the SIR model.Variable nodes are of two types: one type includes pairs of auxiliary variables {(mij, mji)}, the other includes the immunization variables {si}.There is a factor node for each node i of the original graph, including the energetic term Ei and all hard-constraints {Ψ i } involving node i as first label.
discrete variables (i.e.s) defined on the vertices of a graph it is possible to apply the cavity method [47,48] and to develop efficient message-passing algorithms, such as BP and MS.
In the following we present the derivation of the BP and MS equations both for the SIR and SIS models, we discuss the approximation methods employed and the convolution technique used to to solve them efficiently.

A. Derivation for the SIR model
In the SIR model, Equations (5a) define the values assumed by the auxiliary variables m for each configuration of immunized nodes s.The relevant variables on which we must trace to compute the probability weight associated with configuration s are {m ij } (i,j)∈E , with the condition that they satisfy (5a).Let us consider the partition function of the problem where the integration over dm selects (for each s) the unique set of values of m that satisfies the constraints and where the energy is The probability associated with a configuration of immunized nodes s is In order to derive the BP equations, we consider an infinite tree, and we marginalize over all variables but j.In this way we obtain the probability Q j (s j ) for the variable s j , in which the factorization over the neighbors i of j comes from the fact that in a tree each sub-branch is independent of the others, and the partial partition function of the sub-branch is represented by the "cavity marginal" or BP message P ij (m ij , m ji ).The latter denotes the joint probability for the auxiliary variables m ij and m ji in absence of the hard-constraints over node j (i.e.i∈∂j Ψ ji ) and of the energetic term E j ; they satisfy the BP equations Unlike most common cases, the BP message P ij depends on both auxiliary variables defined in the edge (i, j).This happens because both of them enter in the definition of the energetic term (in fact, they are both contributing to m j ).On the contrary the set of hard constraints i∈∂j Ψ ji does not depend on the variable m ij .It follows that in cases in which the energetic term E does not depend on m j (e.g. for = 0 in ( 14)), it can be shown that P ij only depends on m ij .
From the BP marginals we can compute the single-site total probability marginal P i (m i ), that represents the probability distribution of the values assumed by the variable m i on that node.The update rule of the BP equations ( 16) can be represented graphically on a factor graph as shown in Fig. 3.
The cavity message P ij (m ij , m ji ) is a real function defined on the square [0, 1]×[0, 1].Since there is no information on the shape of this function, we proceed discretizing the interval [0, 1] in a number N B of bins, and solving numerically (16) assuming that P ij (m ij , m ji ) is defined on a two-dimensional mesh of N B × N B identical cells.Because of discretization, the hard constraints in (16) are only approximately satisfied.In the following, the integrals over m variables will be replaced by sums whenever we consider discretized versions of the equations.With a naive discretization, it is not easy to keep under control the propagation of the error during the iteration of the BP update rule.We also considered a relaxed version of the constraints for which we can bound the error due to the discretization, however we checked numerically that the two versions of the algorithm give practically the same results.
The computation time of the update rule in the BP equations as written in ( 16) scales as , where k i is the degree of node i, making the trace over all configurations of the neighbors practically unfeasible on most graphs.The computational complexity of the BP update rule can be considerably reduced by exploiting the properties of convolutions of messages.For an arbitrary set D ⊆ ∂i \ j, we define the convolution function where T = k∈∂i [1 − m ki p ki ] so that 0 ≤ T ≤ S ≤ 1.Then, for any disjoint sets D 1 and D 2 , we have One can start from the empty set and add one by one the elements of ∂i \ j, using the convolution rules.Finally, the convolution function over the complete set ∂i \ j can be used to compute the out-coming message as where the two terms, defined as refer to node i being immunized or not (the proportionality symbol as usual means that the message has to be properly normalized).Using the convolution trick, the computational complexity of the update rule on a node of degree k reduces to O(kN 3 B ).The factor N 3 B comes from the computation of the trace over the auxiliary variables {m ki } by means of a convolution function: for each value of the N B values taken by T , we have to sum over all N B values taken by S 1 and by S 2 .
Finally, in order to explore directly the optimal immunization assignments, we can define the MS messages For an arbitrary set D ⊆ ∂i \ j, we define the convolution function where T = k∈∂i [1 − m ki p ki ] so that 0 ≤ T ≤ S ≤ 1.Then, for any disjoint sets D 1 and D 2 , we have The MS equations (in their discretized and efficient version) read where C ij is an (irrelevant) additive constant, and

B. Derivation for the SIS model
In the SIS model, the relevant auxiliary variables by means of which we can connect the configuration of immunized nodes with the corresponding energy value are the mean-field variables {m i } i∈V satisfying (7).As done for the SIR model, we introduce Eqs.(7) as a set of hard constraints in a statistical mechanics model by means of the partition function where now we have and the usual energy term The probability associated to a configuration of immunized nodes s is Factor graph representation of the BP equations for the optimal immunization problem in the SIS model.Variable nodes are of two types: one type includes pairs of auxiliary variables {(mi, mj)}, the other includes the immunization variables {si}.There is a factor node for each node i of the original graph, including the energetic term Ei and the hard-constraint Ψi.
We repeat the same derivation as before, assuming that the graph is an infinite tree.Marginalizing over all variables but j, we compute the probability Q j (s j ) for the variable s j as where the BP messages P ij (m i , m j ) satisfy the BP equations and it represents the joint probability that the mean-field variables on i and j assume values m i and m j in absence of the constraint Ψ j and of the energetic term E j .Figure 4 displays a graphical representation of the BP equations (31) on the corresponding factor graph.The single-site total probability marginal P i (m i ) is given by As already described for the SIR model, the BP equations can be solved numerically using a discretized version of the messages, in which the [0, 1] interval is divided in N B bins.In this way, however, the number of operations required to compute the trace in (31) scales exponentially with the degree of the node, therefore we employ again the convolution method.For an arbitrary set D ⊆ ∂i \ j, we define the quantity Then, for any disjoint sets D 1 and D 2 , we have and the BP equations become Again we can derive MS equations where and C ij is an (irrelevant) additive constant.

V. BP RESULTS ON ENSEMBLES OF RANDOM GRAPHS
On a general graph, the BP equations ( 16) and ( 31) are valid under the hypothesis of fast decay of correlations with the distance or replica symmetric (RS) assumption [47,48].Under this assumption, the statistical properties of the system are described by a unique Gibbs state (i.e.replica symmetry), and the BP equations admit a unique solution.Random graphs are natural benchmark structures for evaluating the quality of the results obtained solving numerically the BP equations with histograms and the performances of the corresponding optimization method.In order to isolate and study the effects of immunization on the statistical properties of epidemic spreading, we consider a completely homogeneous setup: a RRG with uniform values of both spontaneous self-infection and disease transmission along the edges, i.e. q i = q, ∀i ∈ V and p ij = p, ∀(i, j) ∈ E. For the sake of simplicity, we also consider uniform loss parameters and uniform immunization costs, i.e. c i = i = 1, ∀i ∈ V .In the BP approach explained in Section IV, we can give a larger statistical weight to allocations of the immunized nodes that correspond to lower values of the energy.Increasing β for > 0 (see (11)), the distribution Q(s) becomes biased towards immunization sets that generate a smaller expected number of infected nodes compared to random immunizations of the same density.In the limit of β → ∞ the weight is concentrated on the minima of the energy function, i.e. on the optimal immunization sets.In this framework an interesting global observable is the generalization of the quantity f defined in Section III when we perform an average over all possible immunization sets s with the corresponding weight Q(s).We call this quantity f .We can exploit the definition of f in terms of the variables m, and use BP to obtain an estimate of f , where P i (m i ) is given by ( 17) and ( 32) for the SIR and SIS models respectively.The chemical potential µ can be used to control the average fraction of immunized nodes, denoted by v , that is computed from the solution of the BP equations as and for Q i (s i ) we use ( 15) and ( 30) for the SIR and SIS models respectively.It is thus possible to compute f as function of v for a fixed choice of the other parameters.We present here results obtained for infinitely large regular random graphs, obtained using the BP equations in the single-link approximation, i.e. when we solve self-consistently the BP equations assuming all nodes to have essentially the same statistical properties.For both SIR and SIS models, the results of the BP equations in the single-link approximation are shown in Fig. 5 with the choice of parameters q = 0.1 and p = 0.5 and different values of and β.In the case of random immunization ( = 0), the results obtained using BP on infinitely large RRG are compared with the average behavior observed by sampling the solutions of (5a)-(5b) (and ( 7) respectively) and by simulating the SIR (and SIS) stochastic process on finite RRG of N = 10 3 nodes.The latter are obtained sampling over 10 3 configurations of immunized nodes for each value of v .The agreement is very good for both models.
Increasing β for = 1, the solutions of the BP equations show a monotonic decrease in the density of infected nodes.The reduction of the infection level is particularly visible for intermediate values of v , while it becomes almost negligible at small and large density of immunized nodes.The fact that the density v is not directly fixed (as for micro-canonical systems), but it is implicitly varied as an effect of tuning the chemical potential µ and then evaluated from the outcome of the BP equations, is the cause of an undesirable issue at large values of β.It happens that, for large β, the free-energy of the statistical mechanics problem is not convex over the whole interval [0, 1] of values assumed by v .The Legendre transform, implicitly used in the cavity method [48], spontaneously selects the convex envelope of the free-energy, limiting the interval of the possible values assumed by the density of immunized nodes.This is visible in Fig. 5, as increasing β for = 1 the curves get interrupted at some (non-zero) value of v .This means that, for that choice of the parameters, smaller non-zero immunization sets are, in energetic terms, less convenient than no immunization at all (i.e. of the point at v = 0).To eliminate this phenomenon, one can fix the value of v introducing an adaptive external field that is self-consistently adjusted during the iterations of the BP equations [49].Using this technique we were able to explore all values of the density of immunized nodes, although the procedure is not guarantee to converge at all values of β.More precisely, we computed the set of points indicated with black dots in Fig. 5, that correspond to the best results (lowest value of f ) obtained varying β at fixed values of v .We expect these results to be very close to the optimal ones for a large graph.
Upon decreasing the spontaneous self-infection probability q, the effect of the optimization becomes more pronounced.This is shown in Fig. 6, where the results of optimization are compared with random immunization both for q = 0.1 and q = 0.01.In the same Figure, we also plot the entropy S of the immunization sets as function of the density v (computed from the solution of the BP equations as in [48]).The red dashed lines, referring to random immunization, is the standard entropy curve given by S and derived from a binomial distribution of v N immunized nodes among the N possible ones.Red full lines instead represent the entropy curves under strong optimization, i.e. corresponding to low energy states.In this case, the entropy becomes negative for sufficiently small values of v , meaning that no configurations of immunized nodes can be found in the Gibbs state defined at that temperature (this is related to the non-convexity argument explained before).In other words, in that region it is (energetically) more convenient to let the disease spread naturally (with no immune) than  performing immunization.
With the BP analysis we can go beyond the study of the average macroscopic behavior and focus on more informative quantities, such as the distribution P i (m i ).In the absence of immunization, P i (m i ) is a delta function because the equations (5a) (and ( 7)) have a unique solution.After averaging over a set of configurations of immunized nodes, each node i is represented by a distribution P i (m i ) instead of a single value of infection probability m i .In a completely homogeneous setup (i.e.random regular graph and uniform parameters) all nodes have the same statistical properties, therefore we can ignore the node label in the distribution and assume P i (m i ) = P (m), ∀i ∈ V .When the set of immunized nodes is drawn from a uniform distribution (i.e. = 0), we can directly compare the shape of P (m) obtained solving the BP equations with that obtained either by explicitly simulating the stochastic process or by repeatedly solving the mean-field equations with a large sample of configurations of immunized nodes.Examples of the results obtained for the SIR model on RRG of degree K = 4, q = 0.1 and p = 0.5 are displayed in Figure 7.The shape of the distribution P (m) varies considerably for different values of the average immunization density v .The comparison between the results for = 0 and the empirical distributions obtained by sampling solutions of (5a)-(5b) and by means of simulations of the stochastic dynamics is very favorable.The agreement between BP and sampling results seems not to be affected by the discretization of the [0, 1] interval employed in order to solve the BP equations (here N B = 200).Also the agreement with the results of simulations is very good, and it improves increasing the number of configurations of immunized nodes employed in the simulations (here we performed an average over 10 4 realizations at fixed density).
The distribution is always very heterogeneous and the average value (reported in Fig. 5) is not at all representative of the behavior of the system.In general, P (m) is characterized by a series of isolated peaks at low values of m and by a continuous distribution in the bulk.The origin of the delta peaks is strictly related with the local structure of the graph around a node after immunization.For instance, the first peak for m > 0 corresponds to isolated nodes, i.e. nodes that are completely surrounded by immunized ones and thus disconnected from the rest of the graph.Such nodes are infected with probability q, that is exactly the position of the first delta peak.A second peak occurs at m = 0.145, corresponding to the probability of a node being part of an isolated infected dimer.For infected trimers (chain of length three), the central node corresponds to m = 0.187, while the external nodes have m = 0.165.These peaks are visible in small chains, each one corresponding to an isolated peak in the distribution.The continuous bulk of the distribution should be instead identified with a superposition of values due to large-scale clusters of infected nodes.
In Figure 8 we reported a similar plot for the SIS model.The structure of P (m) is generally different, with few isolated peaks in which the weight of the distribution concentrates.This effect could be due to the naive mean-field approximation applied at the level of the equations used as hard-constraints, that is less accurate than that applied in the SIR model.
There is a remarkable difference in the shape of P (m) when a non-uniform weight is associated with immunization sets that generate different levels of infection.In the BP formalism, this is done by increasing from zero. Figure 9a displays the distribution P (m) for the SIR model at a fixed value of v = 0.6, = 1 and different values of β.Increasing β the distribution concentrates on a narrower interval of values of m.The information in these plots is important because it can be used to compute, for a given node in a graph, the probability that such node is infected for a given immunization strategy.One can also ask how many immunization sets exist at a given density v that generate an average infection level of f .For a RRG with degree K = 4 and v = 0.6, this is shown in Fig. 9b, where we plot the entropy of immunization sets of fixed density as function of f .These results are obtained from the solutions of the BP equations increasing β from 1 (for = 1).In Figure 10 we report analogue plots for the SIS model with v = 0.6.For a random allocation of immunized nodes (black circles in Fig. 10a), P (m) exhibits (in addition to a delta in zero) a series of narrow peaks at positive m.Increasing β, those at larger m slowly disappear leaving only a delta peak at m ≈ 0.5.Fig. 10b shows the behavior of the entropy curve at v = 0.6 for the SIS model.In both models, the entropy at v = 0.6 does not vanish continuously when the minimum value of infection f is reached, but it remains considerably large.This is in accord with the behavior of the entropy curves displayed in Fig. 6.A continuously vanishing behavior can be observed instead if we select a density value that falls in the region in which, at large β, the entropy curve becomes negative (see later discussion about MS results).
For q = 0.1, the BP results have been obtained using a discretization of the interval [0, 1] in N B = 200 bins.In general, the smallest possible non-zero value assumed by m is m = q, therefore in order to have a sufficiently good resolution, one should always use N B > 1/q.As the time complexity of the BP and MS algorithms scales as N 3 B , working with histograms when q is small is a problem.A possible solution is that of adopting a mixed representation, such as that exploits the knowledge of the position of the first L delta peaks of the distribution P (m) computed solving the underlying equations for {m i } or m ij on small structures.This approach will be developed in a future work.

VI. COMPARISON BETWEEN DIFFERENT OPTIMIZATION METHODS
A comparison between different optimization methods should take into account both the performances of the optimization and the efficiency of the algorithms.Random regular graphs with uniform parameters are not a particularly good setup to compare different methods because all nodes are equally important.This is visible in Fig. 11a  displays, for the SIR model, the curve f as function of v obtained by several different immunization algorithms: a recalculated degree centrality (green dot-dashed line), recalculated eigenvector centrality (blue line), greedy algorithm (maroon dashed line) and density-constrained simulated annealing (red squares).The results on a RRG with N = 10 3 are compared with the best optimization predicted using density-constrained BP equations on infinitely large RRG (violet crosses).For the same optimization problem, the energy as function of the density of immunized nodes is shown in Fig. 11b (we set µ = = 1.0).Greedy and topologically-based algorithms perform very similarly, and approximately the same result is obtained with fixed-density simulated annealing, even in case of rather slow annealing schedules with up to 10 5 steps between β = 0.1 and β = 10 3 .The BP prediction suggests that (at least for infinitely large graphs) slightly lower energy values could be reached for intermediate values of v .Remarkably, using just 100 bins, the MS algorithm was able to find an immunization set at such lower energies (black vertical cross), that is expected to be the optimal one.We note that the plain MS equations do not always converge.To overcome this difficulty, we employed the reinforcement technique [26,27,50,51].We also ran simulated annealing with > 10 6 steps and without density constraints, and we found the same optimal point.Notice that, although the computational time of MC based methods is comparable with that of MS for graphs of few thousands of nodes, simulating annealing becomes unfeasible for large-scale networks.We also analyzed some epidemic properties associated with the immunization set found by MS and compared them to those expected solving the BP equations in the single-link approximation at the same density v = 0.363.Figure 12a displays the behavior of P (m) computed using BP for different values of and β, and the same quantity obtained solving (5a)-(5b) for the configuration of immunized nodes obtained using the MS algorithm (crosses).The accord between the latter and the BP results for large β values is very good.For the same choice of parameters, Figure 12b shows the entropy as function of the density of infected nodes f (obtained performing BP for increasing values of β).The entropy decreases monotonically and vanishes at f ≈ 0.143, that indicates the minimum level of infection attainable at that density of immunized nodes.The vertical line marking the density of immunized nodes found in the solution obtained using the MS algorithm on a graph of N = 10 3 nodes perfectly matches this lower bound.
We performed the comparison between different optimization algorithms on a real-world network with nonhomogeneous topology.The chosen network, a fully symmetric version of a friendship network among 685 students of a school in the U.S. [52], is of moderate size in order to make possible a direct comparison of the different methods with a simulated annealing with sufficiently slow annealing schedule.Figure 13 shows that, increasing the density of immunized nodes, greedy methods rapidly depart from the optimal curve obtained with density-constrained simulated annealing.On the contrary, heuristic algorithms based of topological metrics performs quite well in finding low energy immunization sets, possibly because of the heterogeneous topological structure of the graph.The optimal point found using the MS algorithm is slightly better than those obtained with all other methods.
When inhomogeneous costs are added, the algorithms based on topological metrics are ruled out.Although they might work very well for some values of µ (the parameter governing the tradeoff between terms in the energy), they do not take into account the contribution to the energy associated with the additional parameters.Energy-based methods are more flexible as they adapt to variations of all parameters appearing in the energy.The greedy algorithm is however limited by the incremental procedure, that finds local minima of the energy with a rather small density of immunized nodes.Adding further nodes to the immunization set, the results tend to worsen.On the contrary, simulated annealing and the MS algorithm are able to find immunization sets with a larger number of immunized nodes but overall less expensive.This phenomenon is observed in the data reported in Figure 14, where we analyzed immunization strategies for the SIR model on the famous Zachary's karate club network of N = 34 individuals.On such a small network, the immunization sets obtained by simulated annealing (performed with 10 5 MC steps between β = 0.1 and β = 10 4 ) are likely to correspond to the global energy minima (for each value of v ).In order to generate a rough energy landscape, we considered immunization costs that are correlated with the degree of the nodes (c i = k i /2, where k i is the degree of i). Figure 14 shows that both for q = 0.1 and q = 0.01, changing the trade-off between the terms of the energy (i.e.increasing µ) the optimal immunization set is found at very different (larger) values of the density v .MS correctly finds the optimal immunization set for all cases under study (black vertical crosses in Fig. 14).

VII. DISCUSSION
In this work, we considered the optimization problem of targeted network immunization against epidemic spreads.The corresponding energy function to be minimized is a tradeoff between the costs of immunizing nodes and the expected extent of the infection.This optimization problem turns out to be computationally worst-case intractable; it was proven NP hard even for very simple stochastic propagation models.For two prototypical models of stochastic epidemic spreading, the SIR and SIS models, we considered mean-field equations that can be used to evaluate the level of infection in the stationary state associated to every configuration of immunized nodes.The original energy function and the problem of finding the optimal immunization set can be (approximately) recast in terms of these mean-field variables in a mixed representation in which the selection of nodes to be immunized is represented by binary variables and the local level of infection (in the stationary state) is expressed by continuous variables in [0, 1].In this formulation, constraints and energy terms are local, allowing the application of the cavity method and the development of efficient message-passing algorithms such as BP.Our results obtained using BP equations on random regular graphs shed light on the statistical properties of immunization sets, uncovering in which regions of the parameter space, and to what extent, targeted immunization is actually more effective than random immunization.The zero-temperature limit of these equations gives the MS algorithm, that can be used to find a solution to the optimization problem.We showed, both on synthetic and real networks, that the MS algorithm outperforms several popular immunization methods based on topological metrics and greedy strategies.The solution found using MS is not guarantee to be optimal, therefore we also performed simulated annealing, that is able to reach the optimum, at least for a sufficiently slow (maybe exponential) annealing schedule.For networks of moderate size, we could compare the lowest-energy immunization set found by MC methods with the solution found by MS.The latter was always at least as good as the former, providing an experimental evidence of the validity of the optimization technique.Moreover, unlike MC-based methods, the MS algorithm scales only linearly with the network size.As a drawback, the algorithm scales as N 3 B , where N B is the number of bins necessary to represent the distribution of real values as histograms.We emphasize that the discretization method used in the current implementation is a very naive and straight-forward one, and there are several ways to considerably reduce N B by adopting more efficient representation of the messages (as explained in Sec.V).For this reason we expect that message-passing algorithms as the one proposed here could be used to study the immunization problem even on very large networks.This will be the scope of future research.
The results of the comparison between different optimization methods on a variety of networks show that, as long as the network and the parameters are sufficiently homogeneous, all methods give approximately the same results.In this case, heuristic strategies based on topological metrics, such as degree-centrality or eigenvector centrality could be preferred, as they are very simple and fast.When the energy function includes inhomogeneous costs, simpler heuristics turned out to be sub-optimal in our experiments.Moreover, greedy-based methods, that take into account the correct energy function, are based on a progressive scheme (as immunized nodes are added incrementally) and often fail to reach the ground state whenever the energy landscape is rugged.This is well demonstrated by results obtained on the small, but representative, Zachary's karate club network.
The method presented here could be applied to a number of other optimization problems including, but not restricted to, other epidemic models in discrete or continuous time, provided that similar fixed-point equations are defined for node or edge variables.The control variables in that case could be a set of node-dependent external parameters.Using message-passing techniques it could be possible to select the configuration of external parameters that corresponds to some desired outcome for the global state of the system.This general formulation opens to a wide spectra of applications in problems involving the control of network dynamics.

FIG. 1 .
FIG.1.Black circles represent the average density f of nodes that have been infected by the final (infinite) time in the SIR dynamics on a random regular graph of N = 10 3 nodes and degree K = 4, as function of the transmission probability p for q = 0.1, 0.01, 0.001 (from left to right).The symmetric bars indicate the fluctuations around the average value computed on 10 4 realizations of the stochastic process.The red full line is computed from the solution of (2)-(3).

4 FIG. 2 .
FIG.2.Black circles represent the average density f of infected nodes in the stationary state (τ = 10 4 steps) of the SIS dynamics on a random regular graph of N = 10 3 nodes and degree K = 4, as function of the transmission probability p for q between 10 −1 and 10 −4 .The symmetric bars indicate the fluctuations around the average value computed on 10 4 realizations of the stochastic process.The red dashed line is computed from the solution of(7).

FIG. 5 .
FIG.5.Results for the analysis of the SIR and SIS models on random regular graphs of degree K = 4, with uniform selfinfection probability q = 0.1 and uniform transmission probability p = 0.5.(a) For the SIR model, we plot the average density of nodes that got infected during the epidemic spreading f as function of the average density v of immunized nodes.BP results on infinitely large graphs are reported for = 0 and β = 1 (black full line) and for = 1 and β = 1 (red dashed line), 5 (green dot-dashed line), 10 (blue double-dot-dashed line), and 20 (violet dot-double-dashed line).Results of sampling over Eqs.(5a)-(5b) (orange circles), corresponding to random immunization, are also displayed.(b) For the SIS model, we plot the average density f of infected nodes in the stationary state as function of the average density v of immunized nodes.BP results on infinitely large graphs are reported for = 0 and β = 1 (black full line) and for = 1 and β = 1 (red dashed line), 5 (green dot-dashed line), 10 (blue double-dot-dashed line), and 20 (violet dot-double-dashed line).Results of sampling over Eqs.(7) (orange circles), corresponding to random immunization, are also displayed.

FIG. 6 .
FIG.6.Results for the analysis of the SIR and SIS models on random regular graphs of degree K = 4, with uniform selfinfection probability q = 0.1, 0.01 and uniform transmission probability p = 0.5.For the SIR model (a-b), we plot the average density of nodes that got infected during the epidemic spreading f as function of the average density v of immunized nodes, both in the random case (open circles) and under optimization (full circles).For the SIS model (c-d), we plot the average density f of infected nodes in the stationary state as function of the average density v of immunized nodes, both in the random case (open circles) and under optimization (full circles).

Figure 7 . 6 FIG. 7 .
FIG.7.Distribution P (m) of the probability m that a node got infected during the SIR process on random regular graphs of degree K = 4 for v = 0.05, 0.2, 0.4, 0.6.Results are computed using BP by means of (17) (black full line) on infinite networks, by sampling solutions of (5a)-(5b) (red dashed lines) and by means of simulations of the stochastic dynamics, both on a finite network of N = 10 3 nodes and with a sample of 10 4 immunization sets.

8 FIG. 8 . 6 FIG. 9 .
FIG.8.Distribution P (m) of the probability m that a node is infected in the stationary state of the SIS process on random regular graphs of degree K = 4 for v = 0.1, 0.4, 0.6, 0.8.Results are computed using BP by means of (32) (black full line) on infinite networks, by sampling solutions of (7) (red dashed lines) and by means of simulations of the stochastic dynamics, both on a finite network of N = 10 3 nodes and with a sample of 10 4 immunization sets.

6 FIG. 10 .
FIG. 10.(a) Distribution P (m) of the probability m that a node is infected in the stationary state of the SIS process on random regular graphs of degree K = 4, q = 0.1, p = 0.5, v = 0.6 and different values of β and .(b) Entropy of immunization sets of density v = 0.6 as function of f (obtained in creasing β for = 1).

FIG. 11 .FIG. 12 .
FIG. 11.Comparison between several immunization methods on a RRG with N = 10 3 nodes and degree K = 4 for the SIR model (with q = 0.1, p = 0.5): (a) average fraction of infected nodes f vs. average fraction of immunized nodes v , (b) energy (for µ = 1.0) vs. average fraction of immunized nodes v .Reported results include: Simulated annealing (red squares), eigenvectors centrality (blue full line), degree centrality (green dot-dashed line), greedy (maroon dashed line) and the best obtained solving fixed-density BP equations (violet crosses) and MS (black vertical cross).Random immunization is also reported (black line for the BP results with = 0 and yellow circles for the results of direct stochastic simulations).

FIG. 13 .
FIG. 13.Comparison between the optimization of immunization sets for the SIR model (with q = 0.1, p = 0.5) on a friendship network of N = 685 individuals: (a) average fraction of infected nodes f vs. average fraction of immunized nodes v , (b) energy (for µ = 0.5) vs. average fraction of immunized nodes v .We report results obtained using: an energy-based greedy algorithm (black full line), fixed-density simulated annealing (red circles), degree centrality (blue dot-dashed line) and eigenvector centrality (green dashed line), MS (black cross).

FIG. 14 .
FIG.14.Energy as function of v for immunization sets obtained using degree (green triangles) and eigenvector centrality (black circles), the greedy algorithm (red squares), and fixed-density simulated annealing (blue crosses), for the SIR model on Zachary's karate club network of N = 34 individuals.The MS results are indicated by black vertical crosses.The immunization cost was assumed to be half of the degree of a node.