FINDING PATIENT ZERO: LEARNING CONTAGION SOURCE WITH GRAPH NEURAL NETWORKS
Year: 2020

Abstract:
Locating the source of an epidemic, or patient zero (P0), can provide critical insights into the infection’s transmission course and allow efficient resource allocation. Existing methods use graph-theoretic centrality measures and expensive message-passing algorithms, requiring knowledge of the underlying dynamics and its parameters. In this paper, we revisit this problem using graph neural networks (GNNs) to learn P0. We establish a theoretical limit for the identification of P0 in a class of epidemic models. We evaluate our method against different epidemic models on both synthetic and a real-world contact network considering a disease with history and characteristics of COVID-19. We observe that GNNs can identify P0 close to the theoretical bound on accuracy, without explicit input of dynamics or its parameters. In addition, GNN is over 100 times faster than classic methods for inference on arbitrary graph topologies. Our theoretical bound also shows that the epidemic is like a ticking clock, emphasizing the importance of early contact-tracing. We find a maximum time after which accurate recovery of the source becomes difficult, regardless of the algorithm used.

heading: 1 INTRODUCTION
text: The ability to quickly identify the origin of an outbreak, or “finding patient zero”, is critically important in the effort to contain an emerging epidemic. The identification of early transmission chains and the reconstruction of the possible paths of diffusion of the virus can be the difference between stopping an outbreak in its infancy and letting an epidemic unfold and affect a large share of a population. Hence, solving this problem would be instrumental in informing and guiding contact tracing efforts carried out by public health authorities, allowing for optimal resource allocation that can maximize the probability of an early containment of the outbreak. Disease spreading is modeled as a contagion process on a network Stroock & Varadhan (2007); Pastor-Satorras et al. (2015) of human-to-human interactions where infected individuals are going to transmit the virus by infecting (with a certain probability) their direct contacts. In general, contagion processes can capture a wide range of phenomena, from rumor propagation on social media to virus spreading over cyber-physical networks Centola & Macy (2007); Baronchelli (2018); Wang et al. (2013); Mishra & Keshri (2013). Therefore, learning the source of a contagion process would also have broader impact on various domains, from detecting sources of fake news to defending malware attacks.
Learning the index case, or patient zero (P0), is a difficult problem. In this paper, we model disease spreading as a contagion process (chains of transmissions) over a graph. The evolution of an outbreak is noisy and highly dependent on the graph structure and disease dynamics. In addition, in real-world epidemics, there is often a delay from the start of the outbreak to when epidemic surveillance and contact tracing starts. Hence, we might only observe the state of the graph at some intermediate times without access to the complete chains of transmission. Furthermore, due to its stochastic nature, the same source node might lead to different epidemic spreading trajectories. Finally, learning P0 from noisy observations of graph snapshots is computationally intractable and the complexity grows exponentially with the size of the graph Shah & Zaman (2011).
Most work in learning the dynamics of a contagion process Rodriguez et al. (2011); Mei & Eisner (2017); Li et al. (2018a) have focused on inferring the forward dynamics of the diffusion. In epidemiology, for example, Pastor-Satorras & Vespignani (2001) have studied learning the temporal dynamics of diseases spreading on mobility networks. The problem of learning the reverse dynamics
and identifying diffusion sources has been largely overlooked due to the aforementioned challenges. Two of the most notable exceptions in the area are “rumor centrality” Shah & Zaman (2011) for contagion processes on trees and Dynamic Message-passing (DMP) on graphs Lokhov et al. (2014) but both require as input the parameters of the spreading dynamics simulations.
Our goal is to provide fresh perspectives on the problem of finding patient zero using graph neural networks (GNNs) Gilmer et al. (2017). First, we conduct a rigorous analysis of learning P0 based on the graph structure and the disease dynamics, allowing us to find conditions for identifying P0 accurately. We test our theoretical results on a set of epidemic simulations on synthetic graphs commonly used in the literature Erdös et al. (1959); Albert & Barabási (2002). We also evaluate our method on a realistic co-location network for the greater Boston area, finding performance similar to the synthetic data. While collecting labeled data to train GNN to find P0 may not be possible, training GNN using simulations on real contact-tracing data can provide a fast method for inferring P0 and help with planning and resource allocation. To the best of our knowledge, our work is the first to tackle the patient zero problem with deep learning and to test the approach on a realistic contact network. In summary, we make the following contributions:
• We find upper bounds on the accuracy of finding patient zero in graphs with cycles, independent of the inference algorithm used. • We show that beyond a certain time scale the inference becomes difficult, highlighting the
importance of swift and early contact-tracing. • We demonstrate the superiority of GNNs over state-of-the-art message passing algorithms
in terms of speed and accuracy. Most importantly, our method is model agnostic and does not require the epidemic parameters to be known. • We validate our theoretical findings using extensive experiments for different epidemic
dynamics and graph structures, including a real-world co-location graph of the COVID-19 outbreak.

heading: 2 RELATED WORK
text: Learning contagion dynamics Learning forward dynamics of contagion processes on a graph is a well studied problem area. For instance, Rodriguez et al. (2011); Du et al. (2013) proposed scalable algorithms to estimate the parameters of the underlying diffusion network, a problem known as network inference. Deep learning has led to novel neural network models that can learn forward dynamics of various processes including neural Hawkes processes Mei & Eisner (2017) and Markov decision processes-based reinforcement learning Li et al. (2018a). Learning forward contagion dynamics have also been intensively studied in epidemiology Pastor-Satorras & Vespignani (2001); Vynnycky & White (2010), social science Matsubara et al. (2012), and cyber-security Prakash et al. (2012). In contrast, research in learning the reverse dynamics of contagion processes is rather scarce. Influence maximization Kempe et al. (2003), for instance, finds a small set of individuals that can effectively spread information in a graph, but only maximizes the number of affected nodes in the infinite time limit. Our problem is more difficult as we care not just about the number of infected nodes, but which nodes were infected.
Finding patient zero In order to find patient zero, we aim to learn the reverse dynamics of contagion processes. Shah & Zaman (2011) were among the first to formalize the problem on trees in the context of modeling rumor spreading in a network. Prakash et al. (2012); Vosoughi et al. (2017) studied similar problems for detecting viruses in computer networks. More recent advances proposed a dynamic message passing algorithm Lokhov et al. (2014) and belief propagation Altarelli et al. (2014) to estimate the epidemic outbreak source. Fairly recently, Fanti & Viswanath (2017) reduced the deanonymization of Bitcoin to the source identification problem in an epidemic and analyzes the dynamics properties. On the theoretical side, Shah & Zaman (2011); Wang et al. (2014) analyzed the quality of the maximum likelihood estimator and rumor centrality, but only for the simple SI model on trees. Antulov-Fantulin et al. (2015) found detectability limits for patient zero in the SIR model using exact analytical methods and Monte Carlo estimators. Khim & Loh (2016); Bubeck et al. (2017) proved that it is possible to construct a confidence set for the predicted diffusion source nodes with a size independent of the number of infected nodes over a regular tree. Our work provides fresh perspectives on the patient zero problem on general graphs based on the recent development of graph neural networks
Graph neural networks Graph neural networks have received considerable attention (see several references in Bronstein et al. (2017); Zhang et al. (2018); Wu et al. (2019); Goyal & Ferrara (2018)). While most research is focused on static graphs, a few have explored dynamic graphs Li et al. (2018b); You et al. (2018); Kipf et al. (2018); Pareja et al. (2019); Trivedi et al. (2019). For example, Kipf et al. (2018) propose a deep graph model to learn both the graph attribute and structure dynamics. They use a recurrent decoder to forecast the node attributes for multiple time steps ahead. Trivedi et al. (2019) take a continuous-time modeling approach where they take the node embedding as the input and model the occurrence of an edge as a point process. Xu et al. (2020) propose a temporal graph attention layer to learn the representations of temporal graphs. However, most research is designed for link prediction tasks and none of these existing studies have studied the problem of learning the source of the dynamics on a graph.

heading: 3 CONTAGION PROCESS AND PATIENT ZERO
text: Finding patient zero means tracing the contagion dynamics back to its initial state and identifying the first nodes that started spreading. Here, we describe the disease dynamics on a network using Susceptible-Infected-Recovered (SIR) and Susceptible-Exposed-Infected-Recovered (SEIR) Kermack & McKendrick (1927) compartmental models that assume that infected individuals develop immunity once they recover from the infections.

heading: 3.1 CONTAGION PROCESSES ON NETWORKS
text: In the SIR model, the population is split into three compartments: susceptible (S) who are susceptible to infection by the disease; infected (I) who have caught the disease and are infectious; removed (R) who are removed from consideration after experiencing the full infectious period.
Continuous time model For a contagion process on a graph G with N nodes, each vertex represents an individual who is in contact only with its neighbors. We can represent the graph using the adjacency matrix A ∈ RN×N , where A[i, j] = 1 if two individuals are connected, 0
otherwise. Let Si, Ii, Ri be the average probabilities of node i being in each of the states, with Si + Ii +Ri = 1. The SIR dynamics on a graph is given by Newman (2018):
dSi dt = −β ∑ j AijIjSi, dRi dt = γIi, dSi dt + dIi dt + dRi dt = 0. (1)
where β is the infection rate per contact and γ the recovery/death rate. The basic reproductive rate of a disease R0 ≡ βλ1/γ is defined as the number of secondary infections created by an index case in a fully susceptible population Keeling & Rohani (2011). Here λ1 is the largest eigenvalue of A. The disease will spread and result in an epidemic if R0 > 1.
Discrete time model We can also use an equivalent discrete time SIR model. Let xti ∈ {S, I,R} be the state of node i at time t. For a susceptible node i, its probability to become infected or removed at time t+ 1 is
P (xt+1i = I|x t i = S) = 1− ∏ j (1− βAijIi(t)) , P (xt+1i = R|x t i = I) = γ. (2)
The SIR model doesn’t account for the incubation period, where an individual is infected but not infectious. This is remedied by introducing an “exposed” (E) state, leading to the SEIR model. For a susceptible node i, the probability to enter the exposed state, and becoming infectious at time t+ 1 is
P (xt+1i = E|x t i = S) = 1− ∏ j (1− βAijIi(t)) , P (xt+1i = I|x t i = E) = α, (3)
An infected node eventually enters the removed state with probability γ, which is the same as SIR (2). (2) and (3) yield (1) for very small β as ∏ j (1− βAijIi) ≈ βAijIi (proof in supp. B).
Finding patient zero Finding P0 can be formulated as a maximum likelihood estimation problem for SIR and SEIR models. Specifically, we observe a snapshot of the state of the nodes at time step t as xt := (xt1, · · · , xtN ), with each node’s state xti ∈ {S,E, I,R}. The problem of finding P0 is to search for a set of nodes Z = {i|x0i = I, i ∈ {1, · · ·N}} which led to the observed state xt. A common approach is to find Z such that the likelihood of observing xt is maximized:
Z? = argmaxZ,|Z|≤kP (xt|Z) (4)
where P (xt|Z) is the probability of observing xt with Z being the P0 set. We assume the number of P0s is no larger than k. Estimating the original state of the dynamic system given the future states requires computing the conditional likelihood P (xt|Z) exactly, which is intractable due to the combinatorics of possible transmission routes.

heading: 3.2 FUNDAMENTAL LIMIT OF FINDING PATIENT ZERO
text: The technical difficulty of finding P0 in SIR and SEIR stems from: (1) presence of cycles in graphs (higher-order transmission) (2) the removed state introducing additional uncertainty about temporal order of infections (3) uncertainty of the exact time step of the observed states. For SI dynamics (i.e. infection is permanent) on trees, existing theoretical results Shah & Zaman (2011); Khim & Loh (2016) have established upper bounds on the detection probability based on an estimator called “rumor centrality”. For graphs with cycles, finding P0 becomes more elusive. We derive the fundamental limit considering the case where at time t = 0 one node, P0, is infected and all of the other nodes are susceptible.
Ambiguity of patient zero on cyclic graphs For graph with cycles, if a cycle is embedded within the infected subgraph, it will reduce the accuracy of predicting P0 because multiple spreading scenarios can lead to the same infection pattern in the cycle. For instance, take a 3-regular tree where the infection has started from the root and spread to some level. If any two branches are connected to make a cycle, both branches become equally likely to have spread the disease downstream. Therefore, having cycles in the infected subgraph reduces the accuracy of finding P0. If O(N) nodes in a graph with cycles are infected, cycles will likely interfere with finding P0. Based on this observation, the following theorem estimates the time horizon beyond which finding P0 becomes difficult.
Theorem 1 (Time Horizon). Assume SIR dynamics (1) on a connected graph of N nodes, starting with a single patient zero. Denoting the adjacency matrix by A and its largest eigenvalue by λ1, the average infection probability, both over nodes and choice of patient zero, 〈I(t)〉 ≡ 〈 ∑ i Ii(t)/N〉P0 becomes O(1) after tmax time steps given by
tmax ∼ logN
γ(R0 − 1) , R0 ≡ βλ1 γ
(5)
Proof: The proof (in Supp B.5.1) follows from the exponential behavior of (1) at small t (Newman, 2018)
Ii(t) ≈ ∑ j exp [t(βA− γI)]ij Ij(0) ≈ exp [(βλ1 − γ)t] ( ψ(1) · I(0) ) ψ (1) i , (6)
and the leading eigenvector being positive ∑ i ψ (1) i ≥ ∥∥ψ(1)∥∥ 2 = 1 (Perron-Frobenius theorem).
The maximum detection accuracy of P0 on a connected graph would decrease with increasing number of cycles in the infected subgraph, GI . In special cases we can estimate the penalty to the accuracy due to cycles. For example, in connected Erdős–Rényi (ER) random graphs Erdös et al. (1959), where each edge has independent an probability p, the number of cycles can be estimated. To get a conservative estimate, we focus on triangles, as they are the most prevalent cycles in ER, to get an upper bound for the accuracy in the following theorem. This ignores cases where the presence of triangles causes downstream error or the error arising from other types of cycles.
Theorem 2 (Detection Accuracy). In contagion process on a connected random graph G, with edge probability p and with infected subgraph GI , the prediction accuracy for P0 is bounded from above
Pmax < 1
3 +
2 3 (1− p)( |GI |p 2 ) (7)
The proof (supp. B.5.2) follows from estimating number of triangles in subgraph GI of a dense ER graph and noting each triangle can drop the accuracy of P0 to 1/3.
Figure 2 shows an example of how this upper bound behaves for different values of R0. The graph is a uniformly connected ER of N = 100, p = 2 logN/N and with γ = 0.4. In conclusion, on graphs with cycles, we expect finding P0 after a time tmax ∼ O(logN) to become difficult. This suggests that to find P0 contact-tracing must be done promptly and in early stages.

heading: 4 FINDING PATIENT ZERO WITH GRAPH NEURAL NETWORKS
text: We propose using GNNs for finding P0 and show that we can improve significantly upon state-ofthe-art methods, e.g. DMP. Moreover, using GNNs gives us the distinct advantage that they are model-agnostic and do not require access to the epidemic dynamics parameters or the time t of the graph snapshot. Our goal is not to propose a novel graph neural network architecture, but to understand the trade-off between different probabilistic inference methods in the context of contagion dynamics. Before we discuss our GNN solution, we briefly review Dynamic Message Passing.
Dynamic Message Passing DMP Lokhov et al. (2014) estimates the probability of every node being the P0 in the SIR model using message-passing equations and approximates the joint likelihood with a mean-field time approach by assuming the following factorization:
P (xt|Z) ≈ ∏
i,xti=S
P (xti|Z) ∏
j,xtj=I
P (xtj |Z) ∏
k,xtk=R
P (xtk|Z) (8)
The algorithmic complexity of the DMP equations over a graph with N nodes and t ≤ T diffusion steps is O(TN2〈k〉) where 〈k〉 is the average degree of the graph. Furthermore, DMP requires providing the SIR epidemic parameters and the time t of the graph snapshot before performing inference. For comparison, on a connected random graph, 〈k〉 > logN , yielding > O(TN2 logN) time complexity for DMP. A trained GNN with l layers (2–8 in our case) has complexity O(N2l) in the inference step and does not require inputting the model parameters. This makes it harder to scale DMP for large or dense graphs. DMP is proven to be exact on trees, e.g. Kanoria et al. (2011), and has been used on more general graphs with reasonable success. In practice, we find that GNN inference can be 50–100 times faster than DMP, as DMP has to exhaustively infer the marginals of S, I,R states for each node for each time t ≥ 0. DMP is trying to accurately infer the distribution over the nodes and auto-regressively infer P0, whereas GNN does not infer node states at each time t and tries to directly infer P0. On the other hand, DMP is an unsupervised inference algorithm and does not “learn”, while GNN is provided with labels in a standard supervised setting.
Relation between Contagion Dynamics and GNNs Our use of GNNs for finding P0 is motivated by the fact that the contagion dynamics (1) are a special case of Reaction-Diffusion (RD) processes on graphs Colizza et al. (2007) which is structurally equivalent to GNNs, as shown next. Proposition 1. Reaction-diffusion dynamics on graphs is structurally equivalent to the messagepassing neural network ansatz.
Denoting pµi (t) ≡ P (xti = µ) of node i being in states such as µ ∈ {S, I,R} or µ ∈ {S,E, I,R} at time t, a Markovian reaction-diffusion dynamics can be written as
pµi (t+ 1) = σ (∑ j F ( Aij · h(pj)µ )) , ha(pi) µ = σ (∑ ν Wµa,νp ν i + b µ ) (9)
where Aaij = θ(Aij)f(A)ij with θ(·) being the step function and σ(·) a nonlinear function. RD on graphs is structurally equivalent to Message-passing Neural Networks (MPNN) Gilmer et al. (2017), as RD involves a message-passing step and a node-wise interaction among features (Supp. B.3),
same as MPNN. We choose the simpler architecture of GCN as in (10) instead of general MPNN. Finding P0 requires learning the backward dynamics of RD, which seems to require the inverse of the propagation rule (PR). Yet, each node can only get infected by its neighbors, so even the backward dynamics requires message passing over the same adjacency matrix and should again have the structure of RD.

heading: 4.1 LEARNING WITH GRAPH NEURAL NETWORKS
text: We employ a state-of-the-art GNN design, suggested by Dwivedi et al. (2020). We make several modifications to the model architecture to fit our problem. Given one-hot encoded node states xti ∈ {0, 1}M as the GNN input, where M is the number of states and where the states are either {S,E, I,R} or {S, I,R}, we first apply a linear transformation h(0)i = Uxti with U ∈ RC×M . Denote the output of layer l by h(l)i , where i is the node index. We use graph convolutional network (GCN) Kipf & Welling (2016) in each layer g(h) = σg (f(A) · h ·W + b), where W ∈ RC×C , b ∈ RC and f is called the propagation rule in GCN. We use f(A) = D−1/2AD−1/2 for the propagation rule, where Dij = δij ∑ k Aik is the degree matrix. To include features of the central node, instead of adding self-loops, we use residual connections between GCN layers and notice a significant increase in model performance. The action of these higher GNN layers is given by
h (l+1) i = h (l) i + σ ( BN(g(h(l)i )) ) , yi = P · ReLU(Q · h(L)i ) (10)
where L is the number of layers and the output layer is parameterized by Q ∈ RD×D and P ∈ R1×D to generate yi ∈ R, representing the probability that node i is P0. BN(·) denotes Batch Normalization and σ(·) is a leaky-relu nonlinear activation function.

heading: 5 EXPERIMENTS
text: We perform extensive studies on performance of our GCNs in finding P0 in SIR and SEIR dynamics over synthetic graphs with various graph topologies. In addition, we generate synthetic epidemic outbreaks on top of a real world co-location network using a SEIR compartmental model that is calibrated to simulate a contagion process with characteristics similar to a COVID-19 outbreak.
Experimental Setup We compare the performance of DMP Lokhov et al. (2014) and different variants of GCNs, following the architecture we described in sec. 4.1:
• DMP: Dynamic Message Passing algorithm Lokhov et al. (2014), we sample a graph snapshot O at time t with each node having a state xti ∈ {S, I,R}, and select the node i that has the highest likelihood of generating O, that is P0= argmaxiP (O|x0i = I). • GCN-S: symmetric GCN Kipf & Welling (2016) f(A) = D1/2AD1/2, GCN-R: random walk f(A) = D−1A, GCN-M: mixture of propagation rules f(A) = A||D1/2AD1/2 • GAT: Graph Attention Network Veličković et al. (2017) We train our models using DGL Wang et al. (2019) with a PyTorch backend. The task is to predict the probability for each node being P0 given a single graph snapshot. We report performance averaged over four random seeds. Details on training and hardware are given in Supp. A.2 We additionally report inference run times.
Evaluation Metrics We use top-1 accuracy to understand the effectiveness of our method. However, due to the ambiguity of detecting patient zero, as elaborated in our theoretical analysis, top-1 accuracy may not be the only evaluation measure to be relied upon. Therefore we also calculate the normalized rank defined byRt = 1− 1 |Dt|N ∑ u∈Dt ru whereDt is the set of test samples at time t, N is the size of the graph and ru is the index of the ground truth P0 in the reverse-sorted probability distribution. Normalized rank is a retrieval metric that tells us how high the correct patient zero was in the learned output distribution. It demonstrates the quality of the output distribution in learning the stochastic dynamics and helps us understand how high was P0 ranked even if it was not ranked the first.

heading: 5.1 EXPERIMENTS WITH SYNTHETIC NETWORKS
text: We first experiment on three synthetic graph models: ER random graph, Barabási-Albert (BA) graph Albert & Barabási (2002) and Random Geometric Graph (Geometric or RGG) Dall & Christensen (2002), and use NDLib (Rossetti et al., 2017) to simulate SIR and SEIR epidemic dynamics on the graph (Supp. A.1). We pick a random P0 seed node at time t = 0 and run S(E)IR a fixed number of steps T . The epidemic parameters (α, β, γ) are chosen such that we can vary R0 to study model performance. We set γ = 0.4 and β = R0γ/λ1 where λ1 is the largest eigenvalue of the graph. For SEIR, we set α = 0.5. We generate 20, 000 simulations and use 80 − 10 − 10 train-validation-test split. For each sample we select t ∈ {1, · · ·T} uniformly at random and try to predict P0 at time t = 0 given the graph adjacency matrix A and node features xti.
We first compare the top-1 prediction accuracy in SIR and SEIR for different models averaged over 1 ≤ t ≤ T . Fig 3 compares the prediction accuracy and average inference time for different models. We can see that GNN-based models outperform the baseline DMP both in accuracy and efficiency. We also want to note that the training time for GNNs is under 7 hours, significantly less than the fastest DMP run of 14.40 hr. It is also worth emphasizing that DMP requires explicit input of β, γ and t while GNNs are model agnostic.
To validate our theory, we plot the theoretical accuracy upperbound and the empirical accuracy obtained from GNN in Fig. 2. We note that the time scale and drop in accuracy are consistent with our theoretical results on tmax (28) and the upper bound on accuracy Pmax (7). Combined with the fact that all our GNN models have comparable accuracies, this suggests that
our GCN-based models may be approaching the fundamental limits we described in 3.2.
Fig. 4 shows the trend of accuracy decay over the time steps t for different graph structures and R0 values. As expected, the accuracy is highest on a tree and when t is small. In graphs with cycles (BA-Dense, ER-Dense, and Geometric) we also observe a nontrivial drop in accuracy which depends both on t and R0. For SIR we observe a drop in accuracy as a function of R0 and time, consistent with our theoretical upper bound. The decay is slower for SEIR as the latent stage adds a delay to the spread of the epidemic. The normalized rank of P0 remains high even over longer time horizons, indicating that P0 could be narrowed down to small subset of the population with impressive accuracy.

heading: 5.2 BOSTON CO-LOCATION NETWORK AND COVID-19 EPIDEMIC TRAJECTORY
text: Our real-world dataset consists of a co-location graph and simulations of an epidemic with the natural progression of COVID-19. The co-location graph is constructed using the Cuebiq data1 for two weeks from 23 March, 2020 to 5 April, 2020 (N = 384, 590 nodes). To reduce computational costs, we sample a subgraph with N = 2, 689 nodes and |E| = 30, 376 edges while maintaining the
1https://www.cuebiq.com/about/data-for-good/ derived from Klein et al. (2020a;b)
degree distribution and connectivity patterns of the original graph. For the epidemic simulations, we run a modified SEIR model with asymptomatic infectious states on the co-location graph with R0 resembling COVID-19 Chinazzi et al. (2020) and accordingly set R0 = 2.5. Each simulation contains 1 patient zero, selected uniformly at random. The simulation is run for 50 days. We create a dataset with 10, 000 samples and an 80− 10− 10 train-validation-test split (supp. C).
The top-k accuracy performance over different days when the graph snapshot was observed are shown in Fig. 5a. We can see that the top-1 accuracy falls steadily over time, the top-(10, 20) accuracy remains fairly high for the first two weeks suggesting that we can retrieve P0 in the most likely 20 nodes out of a total 2, 689 candidates.
Interestingly, while the top-1 accuracy decreases significantly, degrading by 50% after 14 days, using normalized rank, the model can narrow down the set of patient zeros accurately even later in the epidemic, as shown
in Fig 5b. For the normalized rank, P0 can be recovered fairly accurately in the first two weeks of the epidemic. These results highlight an important trade-off between accurately determining patient zero and retrieving the general infected region.

heading: 6 CONCLUSION
text: We study contagion dynamics on a graph using graph neural networks (GNNs) to learn the reverse dynamics of contagion processes and predict patient zero. We evaluate our method against different epidemic models on both synthetic and a real-world contact network with a disease with the natural history and characteristics of COVID-19. We observe that GNNs can efficiently infer the source of an outbreak without explicit input of dynamics parameters. Most notably, GNN accuracy approaches our predicted theoretical upper bound, indicating that further architecture refinements may not improve performance significantly. In addition, GNN is over 100x faster for inference than classic methods for arbitrary graph topologies. Extensions of this work may include learning using sequences of graph snapshots, as well as allowing a set of patient zeros.

heading: A APPENDIX

heading: A.1 DATASET DETAILS
text: We use three graph models: ER random graph, Barabási-Albert (BA) graph Albert & Barabási (2002) and Random Geometric Graph (Geometric or RGG) Dall & Christensen (2002). The density of all three (|E|/ ( N 2 ) ) is adjustable, but BA can produce exact trees. Fixing the number of nodes to N = 1, 000, we first obtain one random instance of tree BA, and dense BA, ER and Geometric graphs with |E| ≈ 10, 000 using the NetworkX library (NetworkX developer team, 2014) and then use NDLib (Rossetti et al., 2017) to simulate SIR and SEIR epidemic dynamics on the graph (supp. A.1). For each sample graph, we pick a P0 seed node i at random to be the patient zero at time t = 0 and then we run S(E)IR a fixed number of steps T . The epidemic parameters (α, β, γ) are chosen such that we can vary R0 to study model performance. We set γ = 0.4 and β = R0γ/λ1 where λ1 is the largest eigenvalue of the graph. For SEIR, we set α = 0.5. We generate 20, 000 simulations and use 80− 10− 10 train-validation-test split. For each sample we select t ∈ {1, · · ·T} uniformly at random and try to predict P0 at time t = 0 given the graph adjacency matrix A and node features xti.
Table 1 describes the details of the synthetic datasets.

heading: Dataset # of Nodes # of Edges Density Diameter

heading: A.2 TRAINING AND HARDWARE
text: We train the model with an ADAM optimizer for 150 epochs with an initial learning rate of 0.003 and decay the learning rate by 0.5 when the validation loss plateaus with a patience of 10 epochs. We perform hyperparameter tuning over a validation set with a random search strategy. We sweep over the hyperparameter space and track our experiments using Weights and Biases Biewald (2020) choosing the model with the lowest validation error. We run our experiments on Nvidia 2080Ti GPUs and report performance averaged over 4 random seeds.

heading: A.3 HYPER-PARAMETER DETAILS

heading: A.4 NOTES ON DMP IMPLEMENTATION
text: We include DMP Lokhov et al. (2014) as a baseline against our proposed GNN based method. As DMP does not have code that is publicly available, we implemented DMP using Python for a fair
comparison with GNNs. Accordingly, our implementation of DMP uses DGL Wang et al. (2019) which enables us to vectorize belief propagation (BP) and marginalization and now it runs in parallel for all nodes.
Given a graph G(V,E), we observe Ot as the state of the graph with nodes i ∈ V . DMP employs MLE estimation to determine the node iP0 that may have led to the observed snapshot O. For a single sample in our dataset D, we use algorithm 1. In order to implement DMP efficiently, we implemented it as a message-passing on a graph using DGL. We sequentially initialize node and edge features for all node i and then as we obtain N = |V | set of graphs with node i acting as P0 in Gi. DMP then allows us to obtain i = argmaxiP (O|i). The advantage of our implementation then is that we can process all N graphs in parallel as if it were one large graph with N2 nodes and E2 edges thanks to DGL’s support for batching graphs. A salient feature of using DGL is that the message passing framework allows us to additionally process all the nodes and edges for a single time step t in parallel. The nature of BP algorithms do not allow us to do away with the for-loop over time t and that remains the only sequential aspect of our implementation. Finally, we use algorithm 1 to process each sample in our test set sequentially. It should be noted that we can further vectorize over a batch of samples in our test set. However, the memory required for DMP is O(bN2E2) with b being the size of the batch and so memory requirements quickly blow up. Accordingly, we leave this aspect of implementation for future work.
Algorithm 1: Dynamic Message Passing given graph G, snapshot O and time t for i ∈ V do
set node i to be P0 initialize node features and edge features with eq (12, 13) in DMP; for (t = 0; t < t; t = t+ 1) do
for e ∈ E do perform message passing with eq (15, 16, 17) in DMP
for j ∈ V do marginalize and update node states with eq (18, 19, 20) in DMP.
Calculate P (O|i) with eq 21 in DMP. return i = argmaxiP (O|i)
A.5 EFFECT OF VARYING NUMBER OF GCN-S LAYERS ON TOP-1 ACCURACY
Fig. 6 shows the top-1 accuracy of P0 of the GCN-S model for varying number of layers. We do not observe a significant effect coming from the number of layers. This may be due to the accuracy limitations with tmax and cycles affecting all the models equally, and superseding other effects such as the diameter of the graph. Another possible reason may be that the 20,000 samples on a graph of 1,000 nodes has many repetitions of the same P0, resulting in both shallow and deep models memorizing patterns.

heading: B THEORETICAL ANALYSIS

heading: B.1 EARLY STAGE EVOLUTION OF SIR AND SEIR
text: The SIR equation on a graph are
dSi dt = −β ∑ j AijIjSi, dRi dt = γIi, dSi dt + dIi dt + dRi dt = 0. (11)
In very early stages, when t 1/γ and ∑ i Ii +Ri N , we have Si ≈ 1 and we have exponential for Ii because
dSi dt = −dIi dt − dRi dt = −dIi dt − γIi
dIi dt = β ∑ j AijIjSi − γIi ≈ ∑ j (βAij − γδij) Ij
Ii(t) ≈ ∑ j (exp[t (βA− γI)])ij Ij(0) (12)
Expanding this using the eigen-decomposition A = ∑ i λiψiψ T i yields eq. (29).

heading: B.2 TRANSITION PROBABILITIES
text: More generally, when the graph is weighted, the probability of susceptible node i getting infected depends on Aij and the probability of node j being in the infected state. For brevity, define pµi (t) ≡ P (xti = µ), with µ ∈ {S, I, ..., R}. The infection probability in SIR (2) can be written as
P (xt+1i = I|x t i = S) = 1− ∏ j ( 1− βAijpIj ) = β ∑ j Aijp I j − β2[ApI ]2 +O(β3). (13)

heading: B.3 REACTION DIFFUSION FORMULATION
text: For brevity, define pµi (t) ≡ P (xti = µ). In a network diffusion process the assumption is that node i can only be directly affected by state of node j if there is a connection between them, i.e. if Aij 6= 0. This restriction means that the general reaction-diffusion process on a graph has the form
Fa(A; p) µ i ≡ ∑ j fa ( ga(A)ijha(pj) µ )
(14)
pµi (t+ 1) = F (A; p(t)) µ i = σ ({Fa(A; p(t)) µ i }) (15)
With
ga(A)ij = θ(Aij)g̃a(A)ij ha(pi) µ = σa (∑ ν Wµa,νp ν i + b µ a ) (16)
where θ(·) is the step function and σa(·) a nonlinear function. In regular diffusion on a graph, we have two states S, I and diffusion is changing the S → I state. The probability Pij ≡ P (xt+1i = I|xtj = S) of node i getting infected at t+ 1, given node j was in the infected state at time t, can be expressed in the form of is determined by the adjacency matrix Aij because node j can only infect its neighbors. The infection probability is given by pIi (t+ 1) = βAijp I j (t) and p S i = 1− pIi . Hence, for diffusion
f1(x) = x g1(A) = βA, h1(pj) µ = ∑ ν δµI δ I νp ν j . (17)
In regular diffusion there is no condition on the target node i and even if it is in the I state the dynamics is the same. In the SI model, however, the infection only spreads to i if it is in the S state. Thus, we have to multiply the dynamics by pSi ≡ P (xti = S) which yields
pIi (t+ 1) = βAijp I j (t)p S i . (18)
This can still be written as (9) by adding the extra functions
f2(x) = x, g2(A) = I, h2(pj) µ = ∑ ν δµSδ S ν p ν j (19)
and having
pi(t+ 1) I = F1(A; p(t)) I iF2(A; p(t)) S i (20)
where Fa = fa(ga · ha) are as in (14). More complex epidemic spreading models such as SIR and SEIR can also be written in a similar fashion. In SIR and SEIR the rest of the dynamic equations are linear and do not involve the the graph adjacency A at all, meaning ga(A) = I in the rest of the equations.

heading: B.4 DISCRETE TIME AGENT-BASED SIR AS A REACTION DIFFUSION SYSTEM
text: The agent-based models (2) and (3), which correct for double-counting of infection from multiple neighbours, are sometimes written as
P (xt+1i = I|x t i = S) = 1− (1− β)ξi , (21)
where ξi is the total number of neighbors j of i which are infected, meaning xtj = I . We will first show that this is a special case of the form given in our paper. First, note that in (2) the terms can also be written as
(1− β)ξi = ∏ j ( 1− βδxtj ,I ) (22)
In the probabilistic model, we have to replace the strict condition of j being in the I state with its probability, so δxtj ,I → P (x t j = I) = p I j (t).
P (xt+1i = I|x t i = S) = 1− ∏ j∈∂i ( 1− βÂijpIj ) (23)
and for small β yield
P (xt+1i = I|x t i = S) = β ∑ j Âijp I j − β2 ∑ j,k Âijp I j Âikp I k +O(β 3) (24)
which yields the simplified equation pi(t + 1)I = pSi (t) ∑ j βÂijp I j (t). Note that if the infection
rate per time step β is large β ∑ j Âijpj can exceed 1, rendering (24) inconsistent with p I i being probabilities. Both (23) and (24) both can be written in the form of RD (15) and (9). We utilize the h1, g1 and h2, g2 found for diffusion (17) and SI (19)
F1(A; p) µ i = ∑ j log ( 1− βÂijh1(pj)µ ) F1(A; p) µ i = h2(pi) µ (25)
and defining the probability as
pIi (t+ 1) = F1 S i
( 1− exp [ F2 I i ]) = pSi (t) 1−∏ j ( 1− βAijpIj (t) ) ≈ βpSi (t)
∑ j Aijp I j (t) (26)

heading: B.5 PROOFS
text: Proposition 2. Reaction-diffusion dynamics on graphs is structurally equivalent of the messagepassing neural network ansatz.
Proof: Analyzing the full stochastic model requires closely tracking the individual events and varies in each run. Hence, we will work with mean-field diffusion dynamics using transition probabilities, instead. Denoting pµi (t) ≡ P (xti = µ) of node i being in states such as µ ∈ {S, I, ..., R} at time t, a Markovian reaction-diffusion dynamics can be written as
pµi (t+ 1) = σ ∑ j F ( Aij · h(pj)µ ) , ha(pi)µ = σ(∑ ν Wµa,νp ν i + b µ ) (27)
where Aaij = θ(Aij)f(A)ij with θ(·) being the step function and σ(·) a nonlinear function. To see this, note that RD processes on graphs involve a message-passing (MP) step (e.g. an infection signal coming from neighbors of a node), and a reaction step where messages of different states µ passed to node i interact with each other on node i. RD dynamics such as the SIR and SEIR models are also Markovian and the probability pµi (t) only depends on the probabilities at t− 1. These are also the conditions satisfied by MPNN. In (27), A are a set of propagation rules for the messages, which are only nonzero where A is nonzero, same as the aggregation rule in MPNN. To have interactions between states µ occurring inside each fixed node i, h(pi) can mix the states µ but not change the node index i, leading to the form of h(pi) in (27), which is the general ansatz for a neural network with weight sharing for nodes, same as in MPNN, and graph neural networks in general.

heading: B.5.1 PROOF OF THEOREM 3
text: Theorem 3 (Time Horizon). Assume SIR dynamics (1) on a connected graph of N nodes, starting with a single patient zero. Denoting the adjacency matrix by A and its largest eigenvalue by λ1, the average infection probability, both over nodes and choice of patient zero, 〈I(t)〉 ≡ 〈 ∑ i Ii(t)/N〉P0 becomes O(1) after tmax time steps given by
tmax ∼ logN
γ(R0 − 1) , R0 ≡ βλ1 γ
(28)
Proof: Consider the spectral expansion A = ∑N a=1 λaψ
(a)ψ(a)T , with λ1 > · · · > λN . In (1) early in the disease spreading, all nodes are susceptible, meaning Si ≈ 1, Ri ≈ 0, and Ii ≈ 1− Si. Thus, combining the three SIR equations, keeping only Ii, the infection spreads as Newman (2018)
Ii(t) ≈ ∑ j exp [t(βA− γI)]ij Ij(0) ≈ exp [(βλ1 − γ)t] ( ψ(1) · I(0) ) ψ (1) i , (29)
Here, I is the identity matrix, λ1 is the largest eigenvalue of A and ψ(1) is the corresponding eigenvector. Averaging over a uniform choice of patient zeros, for the average infection probability we have
〈I(t)〉 ≈ 1 N
exp [(βλ1 − γ)t] 〈 ψ(1) · I(0) 〉 P0 ∑ i ψ (1) i ≥ 1 N exp [(βλ1 − γ)t] (30)
where we used the inequality between L1 and L2 norms to get 〈 ψ(1) · I(0) 〉 P0 = ∥∥ψ(1)∥∥
1 ≥∥∥ψ(1)∥∥
2 = 1. Connectedness means A is irreducible and by the Perron-Frobenius theorem its leading eignvector is positive, hence ∑ i ψ (1) i = ∥∥ψ(1)∥∥ 1 ≥ ∥∥ψ(1)∥∥ 2 . Setting the lower bound of (30) equal to 1 and solving for t we get (28).
B.5.2 PROOF OF THEOREM 2 Ptri
Proof: If P0 is in a triangle, we may miss it 2/3 of the times. Thus, the probability of detecting P0 is bounded by P < 1− Ptri × 2/3, where Ptri is the probability that P0 is in a triangle. Since edges in G are uncorrelated, each having probability p, GI is also a connected random graph with the same edge probability p. Hence, in GI all nodes have degree k ≈ p|GI |. Ptri is one minus the probability that none of the k neighbors of P0 are connected, i.e. Ptri = 1− (1− p)( |GI |p 2 ), which proves the proposition.

heading: C COVID-19 DATA AND SIMULATIONS
text: Geolocation data Mobility data are provided by Cuebiq, a location intelligence and measurement platform. Through its Data for Good program (https://www.cuebiq.com/about/ data-for-good/), Cuebiq provides access to aggregated and privacy-enhanced mobility data for academic research and humanitarian initiatives. These first-party data are collected from users who have opted in to provide access to their GPS location data anonymously, through a GDPR-compliant framework. Additionally, Cuebiq provides an estimate of home and work census areas for each user. In order to preserve privacy, noise is added to these “personal areas”, by upleveling these areas to the Census block group level. This allows for demographic analysis while obfuscating the true home location of anonymous users and preventing misuse of data.
Colocation network The method for constructing the co-location graphs is as follows. First, we split each day into five minute time windows, resulting in 288 time bins per day. For every location event, we use its timestamp to assign it to a time bin, then assign the longitude-latitude coordinate of the observation to an 8-character string known as a geohash. A geohash defines an approximate grid covering the earth, the area of which varies with latitude. The largest dimensions of an 8-character geohash are 38m x 19m, at the equator. If a user does not have an observation for a given time bin, we carry the last observation forward until there is another observation. We finally define two users to be co-located — and therefore to have a timestamped edge in the graph — if they are observed in the same geohash in the same time bin. Accordingly, our co-location graph is constructed by observing the greater Boston area over two weeks from 23 March, 2020 to 5 April, 2020 and results in a graph with N = 384, 590 nodes. To reduce computational costs, we sample a subgraph with N = 2, 689 nodes and |E| = 30, 376 edges with similar degree distribution and connectivity patterns as the original graph and can be observed in Fig 7.
Epidemic simulations in real data. We run a SEIR model on the real co-location network. In doing so, we select parameters and modify the structure of the model to resemble the natural history of COVID-19 Chinazzi et al. (2020). At each time step nodes, according their health status, can be in one of five compartments: S, E, I , Ia, or R. Thus, we split infectious nodes in two categories. Those that are symptomatic (I) and those that are asymptomatic (Ia). The first category infects susceptible node, with probability λ per contact. The second category instead with probability raλ. We set ra = 0.5 and consider that probability of becoming asymptomatic as pa = 0.5. The generation time, that is the sum of incubation (α−1) and infectious period (γ−1), is set to be 6.5 days. Specifically, we fix α−1 = 2.5 and γ−1 = 4 days. In a single, homogeneously mixed, population the basic reproductive number of such epidemic model is R0 = (1− pa + rapa)β/γ where β is the per capita spreading rate Keeling & Rohani (2011). Here however, the epidemic model unfolds on top of the real co-location network. Hence, infected nodes are able to transmit the disease only via contacts (with susceptible individuals) established during the observation period. As mentioned above, the value of R0 is defined by the interplay between the disease’s parameters as well as the structural properties of the network Pastor-Satorras et al. (2015); Masuda & Holme (2017). For simplicity we approximate β = 〈k〉λ, where 〈k〉 is the average number of connections in the network. We obtain λ = 0.073 after solving for R0 and plugging in 〈k〉 = 30376/2689 = 11.29. The simulations start
with an initial infectious seed selected uniformly at random among all nodes. We then read and store the time-aggregated network in memory. The infection dynamics, which are catalysed by the contacts between infectious and nodes, take place on such network. The spontaneous transitions instead (i.e. transition from S toE and the recovery process), take place independently of the connectivity patterns. After the infection and recovery dynamics, we print out the status, with respect to the disease, of each node. Finally, we create a dataset with 10, 000 samples and an 80− 10− 10 train-validation-test split.

Review summaries:
This paper studies the problem of source detection in an epidemics when one observes the underlying graph and a snapshot of the population at a given time i.e. who is infected or not infected. For a SIR (or SEIR) model, the authors propose to use GNN for this task. The learning procedure is then the following: given a fixed graph G, the authors create a dataset of snapshots by running a SIR on G.

Backtracking source of an epidemic (Patient Zero (P0)) is one of the important research topics of the current era that helps efficient resource allocation. Many of the existing works in this domain use graph-theoretic measures or message passing algorithms to tackle this problem. In contrast, this paper uses recently emerged Graph Neural Networks (GNNs) to learn and efficiently locate P0. It models disease spreading as a contagion process over a graph. While considering cyclic graphs, the paper shows upper bound on the accuracy of finding P0 and presents a bound of time horizon after which inference become difficult. Experimental results on different real-world and synthetic networks show interesting results. 

An S(E)IR epidemics propagates on a graph, and the goal is to detect its source (P0) only from the observation of the state (S,E,I,R) of every node of the graph at some time $T > 0$. This version of the source detection problem has been studied first by Shah and Zeman (2011) for SI epidemics, as listed in Section 2. The current paper claims to (i) establish new fundamental limits on this problem, showing in particular that after some time the source detection becomes difficult, and (ii) to demonstrate the ability of graph convolutional networks to solve the problem and validate the results on real data.