We study the problem of estimating the origin of an epidemic outbreak: given a contact network and a snapshot of epidemic spread at a certain time, determine the infection source. This problem is important in different contexts of computer or social networks. Assuming that the epidemic spread follows the usual susceptible-infected-recovered model, we introduce an inference algorithm based on dynamic message-passing equations and we show that it leads to significant improvement of performance compared to existing approaches. Importantly, this algorithm remains efficient in the case where the snapshot sees only a part of the network.  Understanding and controlling the spread of epidemics on networks of contacts is an important task of today's science. It has far-reaching applications in mitigating the results of epidemics caused by infectious diseases, computer viruses, rumor spreading in social media and others. In the present article we address the problem of estimation of the origin of the epidemic outbreak (the so-called patient zero, or infection source -in what follows, these terms are used alternately): given a contact network and a snapshot of epidemic spread at a certain time, determine the infection source. Information about the origin could be extremely useful to reduce or prevent future outbreaks. Whereas the dynamics and the prediction of epidemic spreading in networks have attracted a considerable number of works, for a review see [1] [2] [3] , the problem of estimating the epidemic origin has been mathematically formulated only recently [4] , followed by a burst of research on this practically important problem [5] [6] [7] [8] [9] [10] [11] . In order to make the estimation of the origin of spreading a well defined problem we need to have some knowledge about the spreading mechanism. We shall adopt here the same framework as in existing works, namely we assume that the epidemic spread follows the widely used susceptible-infected-recovered (SIR) model or some of its special cases [12] . The stochastic nature of infection propagation makes the estimation of the epidemic origin intrinsically hard: indeed, different initial conditions can lead to the same configuration at the observation time. Finding an estimator that maximizes the probability of the observed configuration is in general computationally intractable, except in very special cases such as the case where the contact network is a line or a regular tree [4, 6, 11] . The methods that have been studied in the existing works are mostly based on various kinds of graph centrality measures. Examples include the distance centrality or the Jordan center of a graph [4] [5] [6] [7] . The problem was generalized to estimating a set of epidemic origins using spectral methods in [8, 9] . Another line of approach uses more detailed information about the epidemic than just a snapshot at a given time [10] . Note, however, that all the present methods are limited, for instance none of them uses efficiently the information about the nodes to which the epidemic did not spread. In this paper we introduce a new algorithm for the estimation of the origin of an SIR epidemic from the knowledge of the network and the snapshot of some nodes at a certain time. Our algorithm estimates the probability that the observed snapshot resulted from a given patient zero in a way which is crucially different from existing approaches. For every possible origin of the epidemic, we use a fast dynamic message-passing method to estimate the probability that a given node in the network was in the observed state (S, I or R). We then use a mean-field-like approximation to compute the probability of the observed snapshot as a product of the marginal probabilities. We finally rank the possible origins according to that probability. The dynamic message-passing (DMP) algorithm that we use in order to estimate the probability of a given node to be in a given state is interesting in itself. A precursor appeared in [13] in a form averaged over initial conditions, which does not lend itself to algorithmic use. Here we derive and use the DMP on a given network for given initial conditions. If averaged also over the graph ensemble, it can be used to obtain the asymptotically exact dynamic equations of [14, 15] for SIR, or those of [16] for avalanches in the random field Ising model. Note that although DMP bares some similarity with the standard belief propagation (BP) method [17, 18] , it is crucially different from BP in several aspects: it is not derived from a Boltzmann-like probability distribution, and it does not need to be iterated till convergence, instead the iteration time corresponds directly to the real time in the associated SIR dynamics. A nice property that DMP shares with BP is that it is exact if the contact network is a tree. We use it here as an approximation for loopy-but-sparse contact networks in the same way that BP is commonly used with success in equilibrium studies of such networks. We test our algorithm on synthetic spreading data and show that it performs better than existing approaches (except for a special region of parameters where the Jordan center is on average better). The algorithm is very robust, for instance it remains efficient even in the case where the states of only a fraction of nodes in the network is observed. From our tests we also identify a range of parameters for which the estimation of the origin of epidemic spreading is relatively easy, and a region where this problem is hard. Hence, our dataset can also serve as a test-bed for new approaches. Let G ≡ (V, E) be a connected undirected graph containing N nodes defined by the set of vertices V and the set of edges E. The SIR model is defined as follows: Each node i at discrete time t can be in one of three states q i (t): susceptible, q i (t) = S, infected, q i (t) = I, or recovered, q i (t) = R. At each time step, an infected node will recover with probability µ i , and a susceptible node i will become infected with probability 1 − k∈∂i (1 − λ ki δ q k (t),I ), where ∂i is the set of neighbors of node i, and λ ki measures the efficiency of spread from k to i. The recovered nodes never change their state. We assume that the graph G and parameters λ ij , µ i are known (or have been already inferred). The general properties and the phase diagram of this model on random networks were studied in many works, see e.g. [3] and references therein. To define the problem of estimation of the epidemic origin, we consider the case where, at initial time t = 0, only one node is infected (the the 'patient zero', i 0 ), and all others nodes are susceptible. After t 0 > 0 time steps (t 0 is in general unknown), we observe the state of a set of nodes O ⊂ V , and the task is to estimate the location of patient zero based on this snapshot. Let us briefly explain two existing algorithms [4, 6, 7 ] that we will use as benchmarks. The authors of [4, 6, 7] considered only the case when all the nodes were observed, O = V . In appendix B we propose a generalisation of these algorithms to a more general case. The most basic measure for node i to be the epidemic origin is the distance centrality D(i) which we define as D(i) ≡ j∈G d(i, j) δ qj ,I + δ qj ,R /µ j , where the graph G is a connected component of the original graph G containing all infected and recovered nodes and only them, and d(i, j) is the shortest path between node i and node j on the graph G. The ad-hoc factor 1/µ j is introduced to distinguish recovered nodes that for small µ j tend to be closer to the epidemic origin. In the existing works this factor was not present, because [4, 6] treated only the SI model, and [7] considered that susceptible and recovered nodes are indistinguishable. The authors of [4, 6] suggested a "rumor centrality" estimator and showed that, for tree graphs, the rumor centrality and the distance centrality coincide. Another simple but wellperforming estimator, Jordan centrality J(i), was proposed in [7] and corresponds to a node minimizing the maximum distance to other infected and recovered nodes: A node i where J(i) is minimal is known as a 'Jordan center' of G in the graph theory literature. Note that in [7] the Jordan center of only the infected notes was used, hence our implementation uses more information. The core of the algorithm proposed in the present work is DMP, explained in the next section, which provides an estimate of the probabilities P j S (t, i 0 ) (respectively P j I (t, i 0 ), P j R (t, i 0 )) that a node j is in each of the three states S,I,R, at time t, for a given patient zero i 0 . Let us first assume that the time t 0 is known. With the use of Bayes rule, the probability that node i is the patient zero given the observed states is proportional to the joint probability of observed states given the patient zero, P (i|O) ∼ P (O|i). We also define an energy-like function of every node E(i) ≡ − log P (O|i), such that nodes with lower energy are more likely to be the infection source. If one were able to compute P (O|i) exactly, finding i which minimizes E(i) would be an optimal inference scheme of the patient zero. As there is no tractable way to compute exactly the joint probability of the observations, we approximate it using a mean-field-type approach as a product of the marginal probabilities provided by the dynamic message-passing To estimate the value of t 0 , we compute the energy E(i, t) for different possible values t, and choose the value that maximizes the "partition function" Z(t) ≡ i e −E(i,t) . Let us explain the dynamic message-passing equations for the SIR model. The proof that these equations are exact on trees, and the connections to some related existing results [13] [14] [15] [16] [19] [20] [21] [22] , are discussed in appendix A. We first define the message P i→j S (t) as the probability for node i to be in the state S at time t in the cavity graph in which node j has been removed. The quantity θ k→i (t) is the probability for node k not to pass the infection signal to node i up to time t, and φ k→i (t) is the probability for node k to be in the state I and not to pass the infection to node i up to time t. The initial condition are θ k→i (0) = 1, and φ k→i (0) = δ q k (0),I . For more details see 3 appendix A. These messages satisfy the recursion rules: Here ∂i\j means the set of nodes neighboring node i, excluding j. The marginal probabilities that node i is in a given state at time t are then given as Hence the algorithmic complexity for computing the energy E(i) of a given vertex i (and therefore the probability that it is the epidemic origin) is O(t 0 N c), where c is the average degree of the graph. We first test our algorithm on random regular graphs, i.e. random graphs drawn uniformly from the set of graphs where every node has degree c. In all the simulations we consider uniform transmission and recovery probabilities λ ij = λ and µ i = µ. In the first example, inset of Fig. 1 , we plot the energy E(i) resulting from the dynamic message-passing of the nodes for which the probability of being the epidemic origin is finite according to our algorithm, we order the nodes according to the energy value. The true epidemic origin is marked with a red cross. We define the rank of candidates for the epidemic origin to be its position in this ranking (the lowest energy node having rank 0). The main graph of Fig. 1 shows the histogram of normalized ranks (i.e. the rank divided by the total number of nodes that were observed as recovered or infected) of the true epidemic origin as obtained from our DMP inference algorithm, compared to the rankings obtained by distance, rumor and Jordan centralities. The DMP inference algorithm considerably outperforms the three centrality measures, with a comparable computational cost. In Fig. 2 we present the average normalized rank of the true epidemic origin for random regular graphs for the whole range of the transmission probability λ, for different values of the recovery probability µ, and snapshot times t 0 . As an estimation for the spreading time t 0 , we take the one maximizing the "partition function" Z(t) = i e −E(i,t) . The distribution of the estimated time is concentrated at the true spreading time t 0 . We find that for different values of µ, DMP inference outperforms the centrality measures (see, e.g., case (a)), except Inset: An epidemic is generated with recovery probability µ = 1, transmission probability λ = 0.6, a snapshot of all the nodes is taken at time t0 = 8 (in this figure we assume we know the value of t0), 242 nodes are observed to be in the I or R state. The dynamic message-passing is used to compute the energy of every node. This energy is finite for 43 nodes; it is plotted as a function of their rank r. The true patient zero is marked by a red cross, and its rank is r(i0) = 2 in this case. Main figure: an epidemic generated with µ = 1, λ = 0.5, t0 = 5. The histogram (over 1000 random instances) of the normalized rank (i.e. the rank divided by the number of R or I nodes in the snapshot) of the true patient zero is plotted for the dynamic message-passing (DMP) inference, as well as for the distance, rumor and Jordan centrality measures. in a special case (b) (µ = 1, corresponding to the deterministic recovery), in a range of 0.3 < λ < 0.58 where Jordan center is a better estimation. In other cases, however, Jordan centrality is less performant. Note that for µ < 1 Jordan centrality does not distinguish between recovered and infected nodes, which partly explains its rather bad performance in that case. The Fig. 2(c) shows the dependence on the spreading time t 0 for fixed values of λ and µ. Note that DMP remains efficient even for relatively large t 0 , when the centrality algorithms fail to make a prediction. Importantly, in some range of parameters, the average normalized rank of the true epidemic origin is not so close to zero (note that the value 1/2 of the normalized rank corresponds to a random guess of patient zero among all the infected or recovered nodes). The problem of estimating the epidemic origin with a good precision is very hard in these regions. In some cases the information about the epidemic origin was lost during the spreading process. For instance for λ > λ c = µ/(c − 2 + µ) [23] the epidemic percolates at large times t 0 log c N ; then the information about the epidemic origin is lost. On the other hand for t 0 < log c N , the epidemic is confined to a tree network and in this case the inference of the origin is easier, cf. Fig. 2(c) . In Fig. 2, (a) and ( Each data point is averaged over 1000 instances. The plots (a) and (b) represent the dependences of the average rank on the infection rate λ, for the snapshot time t0 = 10 and recovery probability µ: (a) µ = 0.5 and (b) µ = 1. In this figure t0 is inferred by the algorithm. The DMP estimator is compared to the Jordan centrality and the distance centrality estimators. The dotted line shows the average fraction of nodes that were infected or recovered in the snapshot, |G|/N , we use this number to normalize the ranks of the epidemic origin. Plot (c) shows the dependence on the snapshot time t0 for λ = 0.7 and µ = 0.5. The dashed line is the average fraction of nodes that were infected and the dash-dotted line is the average fraction of nodes that were recovered in the snapshot; both are normalized to N . our opinion the systematic comparison presented here is a good test-bed for comparing and improving algorithms. We now show the performance of our algorithm in the case where the snapshot is incomplete: a fraction ξ of nodes is not observed. We compare it to the generalizations of Jordan and distance centralities to this case that we propose in appendix B. Fig. 3 (a) gives the average rank of the true epidemic origin. It shows that, with incomplete snapshots, the DMP inference algorithm outperforms both centralities even in the case where for complete snapshots the Jordan centrality was better. This observed robustness of DMP is a very useful property. In the example depicted in Fig. 3 (b) , we studied the performance of DMP on a non-regular network of the U.S. East-Coast power grid which contains N = 4941 nodes with a maximum degree 19 [24] . We see that the DMP estimator gives better prediction for all range of λ. The results for random scale-free and Erdös-Rényi networks are given in the appendix C. Our algorithm is based on an approximation to the Bayes optimal inference, there are two possible sources of sub-optimality on real networks. The first is the fact that the message-passing equations may lead to errors on loopy graphs. The second is the mean-field-like approximation (1) of the joint probability distribution. We have observed that taking into account the two-point correlation in this approximation does not lead to any improvement in our results. It would be interesting to search for better approximations of the likelihood on a general graph. The approximate solution of dynamics of the SIR model in terms of message-passing equations allowed us to develop an efficient probabilistic algorithm for detecting patient zero. Compared to existing algorithms, it generically (except for a narrow range of parameters) . Recovery probability µ = 1, transmission probability λ = 0.5 and λ = 0.7, only the state of a fraction 1 − ξ of nodes is observed at time t0 = 10, assumed to be known. The rank (averaged over 1000 instances) of the true epidemic origin obtained with our DMP inference algorithm is compared to the distance and Jordan centralities. Right: (b) Normalized rank (averaged over 1000 instances) of the true epidemic origin for epidemic spreading with µ = 0.5, and all nodes observed at time t0 = 10, on the networks of the U.S. East-Coast power grid. DMP inference is significantly better than inference based on distance and Jordan centralities. provides an improved estimate for the source of infectious outbreak. It also performs well when the snapshot sees only a part of the network. One superiority of our approach, compared to previous ones, is that it uses efficiently the information about where the epidemic did not spread. As is usual for Bayes inference approaches, our algorithm is versatile and easily amenable to gener-5 alizations: it can straightforwardly treat the case of dynamically changing networks, the case of multiple sources and also the case where the contact network is not fully known. ACKNOWLEDGMENTS
A dynamic message-passing algorithm is developed to compute a source probability density function in mean-field approximation REF ].
Inferring the origin of an epidemic with a dynamic message-passing algorithm
Phys. Rev. E 90, 012801 (2014)
Physics Computer Science Medicine Biology