Rapid spreading of infectious diseases has had devastating impact on societies throughout human history, but has become more critical in the modern society due to the dense population in urban areas and the increase in human mobility facilitated by the global transportation networks. A recent threat is the spread of the COVID-19 disease caused by the SARS-CoV-2 virus, which has led to a pandemic with severe impact on the public health and the global economy. One prominent feature of this disease is that the pre-symptomatic and asymptomatic viral carriers can spread the disease as well, which poses a big challenge on contact tracing and disease containment [1] [2] [3] [4] . Therefore it is crucial to understand the significance of these undetected transmissions and estimate their impact. Of particular importance are parameter regimes where pre-and asymptomatic infections result in a complete breakdown of our ability to identify infected individuals and contain the spread.
Numerous studies that investigate the spread of the COVID-19 disease aim at predicting the causes of the spreading processes and examine the effectiveness of nonpharmaceutical intervention strategies [5] [6] [7] . It is common to model the evolution of the population mass of each group (e.g. susceptible, exposed, infected and recovered) by deterministic differential equations [4, 5, 8] . While being simplistic and tractable, such method assumes homogeneous mixing of the population (in a city or within an age group) and neglects the social contact network structures of the specific instance investigated [9] . Large scale agent-based simulations are also widely used, which provide a more detailed picture of the spreading processes but are very computationally demanding and * b.li10@aston.ac.uk † d.saad@aston.ac.uk suffer from a lack of a principled understanding [10] [11] [12] [13] .
To obtain a reliable statistical description of the process one has to increase the number of samples significantly as the system size increases, which makes the computation prohibitive for large systems. Analytical treatments to the epidemic spreading processes on heterogeneous contact networks are valuable both in providing solutions in specific instances and in exploring the typical macroscopic behavior of ensemble of systems with similar characteristics; the latter also results in generic and intuitive understanding.
In this work, we analyze diseases spreading with presymptomatic transmission such as COVID-19 by studying a variant of stochastic susceptible-exposed-infectedrecovered (SEIR) model on contact networks, in which nodes in exposed state can also spread the disease without showing symptoms. For simplicity, the contact networks are viewed as static, serving as substrates on which the disease spreads. We derive the dynamic messagepassing equations for this model, which provide a good approximation to the complex stochastic spreading dynamics on general networks and facilitates theoretical analyses [14] [15] [16] [17] . Based on this framework, we derive the epidemic thresholds and their dependence on different intervention methods. The emphasis of this work is not on predicting the exact dynamics of the COVID-19 epidemic, but to obtain a more intuitive physical pictures and clarify the effects of some containment strategies.
The remainder of the paper is organized as follows. We introduce the SEIR model in Sec. II, and derive the dynamical equations in Sec. III. We then perform a linear stability analysis of the dynamical equations in Sec. IV, based on which the epidemic thresholds are obtained and analyzed in Sec. V. In Sec. VI, we show that non-backtracking centrality can be used to predict the outbreak profile. In Sec. VII, we investigate the effects of reducing contacts on slowing down the spread of the disease. Finally, we summarize our findings and discuss some limitations and outlook.
The contact network is represented by a graph G(V, E) where V is the set of nodes and E is the set of edges. We assume that the network has only one connected component. Each individual resides on a node, assuming one of four states, susceptible (S), exposed (E), infected (I) and recovered (R) at any particular time step. We assume that a node in the exposed state has contracted the disease but has not developed symptoms yet. Unlike the usual SEIR model [18] , the exposed nodes can also spread the disease. The dynamical process of the modified SEIR model in discrete time is defined in the form of transition probabilities of states of neighboring nodes (say i and j) and the state evolution of an individual node i
where α ji (β ji ) is the probability that node j being in the exposed(infected) state transmits the disease to its susceptible neighboring node i at a certain time step. We assume that each time step corresponds to one day, keeping in mind that a finer time scale can also be considered. At each time step, an existing exposed node i becomes infected (i.e., develops symptoms) with probability ν i , while an existing infected node i recovers with probability µ i . Therefore, the average periods of incubation and recovery are 1/ν i and 1/µ i , respectively. At a certain time step, the symptom-development and recovery processes are assumed to occur after possible transmission activities. Since we will contrast the properties of the SEIR and SIR models, we also introduce the transition probabilities of the latter S(i) + I(j) βji − − → I(i) + I(j),
which has been widely studied in the literature [9] . We remark that both models are Markovian processes, implying an exponential distribution for both symptomdevelopment and recovery times, which may not be fully realistic for many diseases including COVID-19 [9, 19] . Nevertheless, they both represent relevant models that are amenable to analysis and would provide important generic insight. The epidemiological parameters depend on the nature of the disease and the intervention strategies being imposed; they are usually estimated based on observations and can be subject to a high degree of uncertainty. As for COVID-19, the average incubation period is about 5.2 days [3, 20] . Infectiousness is estimated to start from 2.3 days before the onset of symptoms [3] , while it is argued [21] that infectiousness can start much earlier (one needs to look back at 5 days to catch 97% of pre-symptomatic infections). The time needed for the symptoms to disappear depends on disease severity of the patient. In [3] , it is inferred that infectiousness declines rapidly within 7 days. In [22] , it is found (from patients with mostly mild COVID-19) that viral subgenomic RNA, which provide evidence of replicative intermediates of the virus, were detectable up to 8 days after the onset of symptoms. In this work, we define the recovered (R) state where the exposed/infected individual effectively looses infectiousness, irrespective of whether the symptoms persist or not. Therefore, an infected individual who has been put into isolation and can no longer infect others is also categorized to be in state R. To address the COVID-19 disease, we set ν i = 1/5, µ i = 1/8 according to these previous findings.
The transmission probabilities α ji , β ji are more difficult to estimate. In [2] , the transmission rate of undocumented infections was estimated to be 55% of that of documented infections. Based on this finding, we set α ji = β ji /2. For simplicity, we also assume that the parameters are homogeneous, i.e., α ij = α, β ij = β, ν i = ν, µ i = µ. Our framework can accommodate any infection/recovery parameter distributions, but the insight gained from such detailed parameter sets are more limited. Various intervention strategies have different impacts on these epidemiological parameters. Some studies report cases where infected individuals remain asymptomatic throughout the course of the infection; however, these cases are the minority and seem to have a much lower secondary attack rate [23, 24] . To keep the analysis simple, we do not consider asymptomatic individuals who do not become infected prior to recovery but briefly discuss how the frameworks used could accommodate such cases in Appendix B.
Since the exact solutions of the stochastic spreading processes Eq. (1)-(4) are difficult to obtain, various approximation methods have been developed to tackle such complex dynamics [9] . A simple method is the individualbased mean-field (IBMF) approach [9, 25] , which expresses the evolution of the marginal distribution P i σ (t) that each node i belongs to state σ by assuming the independence on the probabilities of neighboring nodes. Consider the SEIR model defined in Sec. II. For node i being in the susceptible state S, it will remain in state S in the next time step if none of its neighbors transmits an infection signal to node i. In the IBMF framework, this occurs with probability k∈∂i 1 − α ki P k E (t) − β ki P k I (t) , where ∂i denotes the set of nodes adjacent to node i. Therefore the evolution of the marginal probability node i in state S is given by
The probability of node i in the exposed state E increases if there is an infection signal from its neighbors, while it decreases with rate ν i (probability of transforming into state I) as infection symptoms appear. The corresponding IBMF dynamical equation is
Similarly, the evolution of P i I (t) and P i R (t) are given by
This approach has been used to investigate similar models addressing the COVID-19 pandemic [26, 27] . However, the drastic simplification based on the independence assumption of probabilities may lead to large approximation errors [25] . One source of errors comes from the mutual infection effect due to this decorrelation assumption [17, 28] . For instance, suppose that a node i, having probability P i E (t) in the exposed state, infects its susceptible neighboring node k at time t, then node k can also reinfect node i at time t + 1 with some probability, which is an artifact of neglecting the correlation between nodes i and k. Such effects need to be correctly accounted for in order to improve accuracy.
The dynamic message-passing (DMP) approach, an algorithm that originates from the statistical physics literature [14] [15] [16] avoids the mutual infection effect by considering the irreversible complete trajectories of the system. Formally, the DMP equations can be derived from the belief propagation equations of dynamical trajectories, which is especially useful when the correct set of dynamical variables are difficult to determine straightforwardly [16, 29] . In this section, we provide an intuitive derivation of the DMP equations of the SEIR model, while we give the more formal derivation based on belief propagation in Appendix A.
Similar to the IBMF approach, the DMP method aims at deriving the evolution of the marginal distributions. Consider the marginal probability P i S (t) of node i being at state S at time t, it is given by
where θ k→i (t) is the probability that node i has not contracted the disease from node k up to time t. We have made the assumption that the probability which node i has not contracted the disease from its neighbors up to time t factorizes as k∈∂i θ k→i (t). This assumption is valid in tree graphs but constitutes a good approximation in many loopy networks [16] .
The message θ k→i decreases if node k transmits the infection signal to node i, which occurs with probability α ki if node k is in state E or with probability β ki if node k is in state I. Therefore, it follows the update rule (12) where ψ k→i (t) is the probability that k is in state E but has not transmitted the infection signal to node i, and φ k→i (t) is the cavity probability (on a graph where node i is absent -a cavity) that k is in state I but has not transmitted the infection signal to node i up to time t.
The message φ k→i decreases if node k transmits the infection signal to node i or changes from state I into state R; note that the two processes can occur at the same time step. On the other hand, it increases if node k changes from state E into I. Therefore, it is updated according to
Similarly, the message ψ k→i (t) decreases if node k transmits the infection signal to node i or changes from state E into state I, while it increases if node k changes from state S into E. In computing the probability increment due to the latter case, one needs to exclude the contribution from node i to node k in the previous time steps in order to avoid the effect of mutual infection. This is achieved through defining
which is computed in the cavity graph assuming that node i has been removed. Then the message ψ k→i (t) follows the update rule
Upon computing the messages {θ k→i (t), φ k→i (t), ψ k→i (t), P i→j S (t)} from the update rules, the marginal probability P i S (t) can be obtained by Eq. (11) . The marginal probabilities of other states can also be determined as
We assume that the nodes are either in susceptible or exposed states at time t = 0, such that the initial conditions are solely determined by P i S (0) as
If node i is selected as the initial exposed node, the initial condition is simply set as P i
Based on the initial data, the messages are solved by updating the DMP Equations (12)- (15) forward in time, after which the marginal probabilities are determined by Equations (11) , (16) , (17) and (18) . The computational complexity of the DMP algorithm is linear in number of time steps and in number of edges as O(|E|T ). Therefore, the DMP approach saves a significant amount of computational resources compared to Monte Carlo (MC) simulations, which requires many realizations to obtain reliable results. On the other hand, it is more demanding than the IBMF approach which deals with node-wise variables rather than edge-wise variables. Nevertheless, if the network is sparse, i.e., the average degree d = 2|E|/N ≪ N , then the DMP approach has the same order of computational complexity as the IBMF approach. This is relevant to the case of disease spreading as the number of close contacts each person has is limited [30] , except for super-spreaders.
Here we evaluate the effectiveness of the developed theories on contact networks, which are either artificially generated or adapted from some realistic human contact data. The realistic contact networks are taken from data sets obtained in the SocioPatterns collaboration [31], where the temporal face-to-face human contacts are projected to static contact networks as described in Appendix C.
Many realistic contact networks exhibit community structures. For instance, in workplaces, people usually interacts more frequently with other people from the same department compared to those from other departments; in schools, students from the same class also interacts more frequently. An exemplar contact network in the workplace (WP2015) from the SocioPatterns data is depicted in Fig. 1(a) . We run the DMP algorithm for this contact network by randomly selecting 5 nodes as the initial exposed nodes and iterating for T = 100 time steps. The evolution of the population of each compartment is shown in Fig. 1(b) . The number of exposed or infected cases first rises and then decreases, and eventually dies out when herd immunity is reached. The trajectories predicted by the DMP approach match well with those from MC simulation in this case. , there is only one exposed node. The accuracy of the DMP approximation is much better than the IBMF approach.
As another example, we evaluate the efficacy of our framework on random regular graphs, where all nodes have the same degree and are connected randomly. Only one node is selected as the initial exposed node, and the system is simulated for T = 30 time steps. The results are shown in Fig. 2 , exhibiting a much better approximation accuracy of the DMP approach compared to IBMF.
In Fig. 3 , we systematically compare the results between theories and MC simulation, showing that DMP provides a much better approximation than the IBMF approach. It is found that the prediction errors of both theoretical approaches generally increase with the epidemiological parameter β, and also depend on the number of initial exposed nodes. Intuitively, large values of β Figure 3 . Comparison between theory and MC simulation for different parameters, where the prediction error is calculated as σ∈{S,E,I} i,t |P i MC,σ (t) − P i theory,σ (t)|. The parameters used are T = 30, ν = 0.2, µ = 0.125, α = β/2. Each data point is average over 5 instances with different sets of randomly selected initial exposed nodes. (a) Random regular graph with N = 100, d = 10. (b) Contact network in the workplace taken from the SocioPatterns data (WP2015). The prediction errors depend on the number of initial exposed nodes and the epidemiological parameter β.
lead to larger growth rate of the infections, in which case a small approximation error in early time steps could be amplified in late times. For large β, the prediction errors typically decrease as the number of initial exposed nodes increases. One possible reason is that the infections spread out from a unique source are correlated, so that the independence assumption in the IBMF and DMP approaches deteriorates [32] , e.g., the assumption that the messages in Eq. (11) factorize does not hold strictly. On the other hand, if there are multiple initial seeds that trigger the outbreak, the infection signals into a node will be less correlated, which partly preserves the decorrelation assumption.
The network topology also impacts on the approximation accuracy of the theories as shown in Appendix D. As mentioned before, DMP avoids the mutual infection by excluding one step backtracking interaction, it does not take into account mutual infections due to backtracking of multiple steps, which can be non-negligible in networks with many short loops. In Appendix D, we find that the approximation accuracy deteriorates significantly when the localization of non-backtracking centrality is present, an effect we will discuss below. Due to this effect, both the IBMF and DMP tend to overestimate the outbreak.
The fate of the spreading processes depends on the epidemiological parameters. For large transmission probabilities β and α, the disease can spread out to a large fraction of the network, while it tends to die out quickly for small transmission probabilities. There exist thresholds for these parameters, above which the epidemic outbreaks occur. One commonly used method to determine epidemic thresholds is to examine whether the disease-free state is linearly stable to small perturbations [9, 33] .
Specifically, the initial disease-free state is perturbed infinitesimally as P i S (0) = 1 − ǫ i ; if such perturbation diverges, then the outbreak tend to spread out globally. At the initial stage, the message θ k→i (t), which denotes the probability node k has not passed the infection signal to node i, can also be expressed as θ k→i (t) = 1 − δ k→i (t). At time t + 1, we have
which implies
where δ k→i (t), ψ k→i (t) and φ k→i (t) have small values. Expanding Eq. (14) and keeping first order of ψ k→i (t) and φ k→i (t) leads to
Then ψ k→i (t + 1) in Eq. (15) can be approximated as
Equations (26) and (13) constitute a linear dynamical system of the messages {φ i→j (t), ψ i→j (t)}.
In the following, we use homogeneous parameters α ij = α, β ij = β, ν i = ν, µ i = µ. To make the linearized dynamical equations more compact, we introduce the 2|E| × 2|E| non-backtracking (NB) matrix with elements
which are non-zero if and only if the directed edge i → j follows right after edge k → l, i.e., in a configuration like k → l(= i) → j, but edge i → j does not backtrack to node k [34] . Then Eqs. (26) and (13) can be written in the matrix form as
where the J is the Jacobian matrix of the dynamical system defined as
where I is the 2|E|-dimensional identity matrix. The spectral radius ρ(J ) of the Jacobian J determines the growth rate of the fastest mode of the linearized dynamics. The appearance of the NB matrix B in the linearized dynamical equation is rooted in the fact the one-step backtracking infection is excluded in the DMP equations (e.g., through Eq. (14)), therefore it also appears in other algorithms for complex networks based on linearizing belief propagation, such as the applications in community detection [35] and percolation [36] .
The condition that ρ(J ) ≥ 1 corresponds to an exponential growth of the linearized dynamics in Eq. (28) , which implies that the disease is likely to spread out globally. The solution of ρ(J ) = 1 marks the phase boundary of the epidemiological parameters. Since the matrix elements of J are non-negative, the Perron-Frobenius theorem asserts that (i) its leading eigenvalue λ max J (defined as the eigenvalue having the largest real part) equals its spectral radius ρ(J ) and therefore is real and non-negative; (ii) there exists an eigenvector with non-negative and nonzero elements corresponding to λ max J [37] . Consider the eigenvalue equation of the Jacobian matrix J
which can be simplifies to
is an eigenvalue of the NB matrix B with eigenvector u, denoted as λ B . The Perron-Frobenius theorem also guarantees that the leading eigenvalue λ max B of B is real and non-negative. It can be shown that the leading eigenvalue
In this way, we relate the dynamical properties of the SEIR model to the epidemiological parameters and network structure properties, where the latter is subtly conveyed through the eigenvalue λ max B . This is in contrast to the growth rate given by the commonly used basic reproduction number R 0 , defined as the expected number of secondary cases caused by a single randomly-selected exposed individual when the rest of the population are susceptible. In the SEIR model considered here, the R 0 is estimated to be (see Appendix F)
which only depends on the averaged degree d , but neglects possible higher order structures of the contact networks. It has long been recognized that the R 0 measure is deficient in network epidemiology [9, 25] . The computation of the leading eigenvalue λ max B can be demanding as the NB matrix is of size 2|E| × 2|E|, which is a large matrix if the underlying network is relatively dense. It has been observed that the spectrum of B can be obtained from a much smaller matrix of size 2N × 2N [35, 38] 
where I N is the N -dimensional identity matrix, D is the diagonal matrix of node degrees with elements D ij = d i δ ij , and A is the adjacency matrix with elements sat-
Intuitively speaking, the reduction of complexity comes from compressing the edge-based data, e.g., φ i→j (t) and ψ i→j (t), to node-based data [35] as shown in Appendix G. This allows us to work with networks of relatively large sizes.
Equation (33) gives rise to the epidemic threshold as predicted by the DMP approach through the solution of λ max J (β, α, ν, µ, λ max B ) = 1, keeping in mind that λ max
The same derivation can also be applied to the IBMF equations, as shown in Appendix G, leading to
, which relates the maximal eigenvalue λ max J MF of the Jacobian matrix J MF of the IBMF equations in Sec. III A and the adjacency matrix A of the network. Solving λ max J MF (β, α, ν, µ, λ max A ) = 1 yields the epidemic threshold as predicted by the IBMF approach.
The epidemic thresholds obtained by the two theoretical approaches are to be contrasted with those obtained from numerical simulations, where the large time limit is taken such that the outbreaks saturate and the final state of each node i is either susceptible or recovered. The fraction of nodes that have been infected is an order parameter r = i P i R (∞)/N defining the phase transition from localized infections to global epidemics. Since statistical fluctuation is large near criticality, one can estimate the critical point through the variability measure [39] 
which peaks at the critical point.
As an example, we consider random regular graphs of degree d = 10, where the leading eigenvalues of the matrices A and B have exact expressions as λ max
irrespective of the network size, as described in Appendix G. We also fix the values of ν and µ, let α = β/2 and consider the phase transition by varying β. It is shown in Fig. 4 (a) that a significant fraction of the systems nodes are affected by the epidemic outbreak above the critical point β c . In Fig. 4 (b), we pinpoint the critical point β c from MC simulations through the variability measure C r , and compare them to those obtained via the IBMF and the DMP approaches, where it is observed that the DMP approach provides a much better estimation. In Appendix G we observed that the epidemic has a small probability to die out even for large β, a behavior that also appears in the SIR model [39] and may impact on the estimation of β c through simulations [28] . This is not captured by the theories which only consider averaged quantities. It would be an interesting future direction to study the deviations from the mean behaviors [40] and possible heterogeneous structures [41] .
We examine the phase diagrams in parameter subspace spanned by β and α in different scenarios via the DMP approach in Fig. 5 . We expect that various social distancing measures effectively reduce β and α, though it is difficult to quantify the exact reduction due to a particular policy. The critical line separating the parameter regions of localized infections and global outbreaks, obtained by solving λ max J (β, α, ν, µ, λ max B ) = 1 for α, has the following expression 
The systems are initiated with 5 exposed nodes and are simulated for a sufficiently long time such that every node that has contracted the disease has recovered. In both (a) and (b), the three vertical lines correspond to the critical points βc obtained via different approaches, and different lines correspond to networks of different sizes. In both the IBMF and the DMP approaches, βc is the same for networks of different sizes as the degree is fixed. In MC simulation, βc is obtained through the largest network (N = 6400) considered. (a) Order parameter r = i P i R (∞)/N as a function of β, obtained in MC simulations. (b) Variability of the order parameter, i.e., Cr = r 2 − r 2 / r , as a function of β. Since the variability Cr obtained in simulations (data points in circle markers) is subject to fluctuations, we apply curve fitting to the experimental data using Gaussian processes with a radial basis function kernel (the fitted curves are shown in solid lines) [42] . The maximum point of Cr in the fitted curve (diamond-shape marker) is marked as the critical point obtained from MC simulations. Inset of (b) shows that the critical point βc obtained in MC simulations approaches the one obtained by the DMP approach as the network size increases.
While there is no pre-symptomatic transmission in the SIR model defined by Eqs. (5) and (6), we can compare with the expression obtained for the critical transmission probability [28] β SIR
As a comparison, we sketch the phase boundaries α c (β) of the SEIR model and β = β SIR c of the SIR model for certain ν, µ and λ max B in Fig. 5(a) . The epidemic will not spread globally when the parameters are in both regions (I) and (II) in the SIR model. In contrast, the disease will die out only in region (I) in the SEIR model. Particular caution is needed in region (II), where β is small enough so that it is safe for the SIR model, but the presymptomatic transmission is sufficiently significant (i.e., α > α c (β)) to cause a global epidemic.
Consider the special case α = 0, solving α c (β) = 0 gives
which coincides with β SIR c . It indicates that the intersection of the two critical lines occurs at α = 0 (as shown in Fig. 5(a) ) is a general phenomenon. Physically, this special case corresponds to the traditional SEIR model where the exposed nodes are not infectious, which has the same epidemic thresholds with the SIR model, irrespective of the incubation period 1/ν.
Similarly, consider the special case β = 0, then
which is effectively the critical point of an SIR model viewing the state I as "recovered", as expected. Equations (40) and (41) explain some of the behaviors in Fig. 5 (b)(c)(d), as explained below. When ν and λ max B are fixed, varying µ will only affect the intersects of the phase boundary on the β-axis (β SEIR c | α=0 depends on µ) but not the α-axis (α SEIR c | β=0 is independent of µ), as seen in Fig. 5(b) . Physically, since different values of µ represent different recovery rates, increasing µ can be effectively realized by isolating all nodes in state I before they recover. If such a policy can be executed strictly and timely, it corresponds to a large µ that can significantly expand the parameter region of the epidemic-free phase, leaving a lot of flexibility in implementing social distancing measures (i.e., many usual social interactions can still be allowed). Nevertheless, if the pre-symptomatic transmission probability α is large enough, e.g., α > α SEIR c | β=0 , then isolating the patients in state I alone is insufficient to slow down the spread; in this case, identifying nodes in state E through contact tracing or mass testing, and/or implementing stricter social distancing measure become necessary.
Similarly, when µ and λ max B are fixed, varying ν will only impact on the intersects of the phase boundaries with the α-axis but not those with the β-axis, as seen in Fig. 5(d) . Physically, different ν values correspond to diseases with different incubation periods. For smaller ν, the exposed nodes have a longer time to infect their neighbors, which makes it more difficult to combat the epidemic spreading.
Finally, although the phase boundary α c (β) of the SEIR model is in general nonlinear, the cases considered in Fig. 5 exhibit an almost linear relation, except for a very small λ max B . This can be seen more explicitly in the limit of large λ max B , where one of the sufficient conditions is that the network has a large average degree (e.g., λ max B = d − 1 in random regular graphs), and requires the transmission probability β and α to be small enough for the disease to die out. Under the condition of large λ max B and small β, Eq. (38) can be approximated as
Equation (42) reducing λ max B can be achieved by limiting the number of contacts between nodes, as will be shown in Sec. VII.
A similar linear relation of the phase boundary can be obtained by the condition R 0 = 1, where R 0 satisfies Eq. (34), yielding
which coincides with Eq. (42) if we identify λ max B as d . It turns out that the condition λ max B = d holds approximately for Poisson random graphs [43] , as will also be shown in section VII. Therefore, the critical line obtained via the basic reproduction number is a good estimation in dense Poisson random graphs, but may become a poor approximation otherwise.
Similar to the eigenvalues, the eigenvectors of the Jacobian matrix also provide valuable information on the dynamics. Consider the eigen-decomposition of the Jacobian as J = a λ a J ξ a (ξ −1 ) a , and note that the messages can be decomposed using the eigenvectors as bases, i.e., ψ(t), φ(t) = a c a (t)ξ a . In light of the linearized dynamics, the component c max (t) corresponding to the leading eigenvalue λ max J , will dominate when λ max J > 1 as the system evolves. Therefore, we can use the leading eigenvector ξ max of the Jacobian J to predict the outcome ψ(T ), φ(T ) of the dynamics.
According to Eq. (31), the two components of ξ max = (ũ,ṽ) are proportional to each other, i.e.,ṽ ∝ũ. Thus, it is sufficient to examine one component only. In what follows, we considerũ, which is the leading eigenvector of the NB matrix B according to Eqs. (32) and (33) . Furthermore, we are ultimately interested in the marginal probabilities, which relate to the incoming messages to each node as seen in Eq. (11) . In the linearized dynamics, the probability of newly infection of node i is given by
Therefore, we need to consider the incoming vector of the leading eigenmodeũũ
which is known as the non-backtracking centrality [43] .
Interestingly, in addition to the leading eigenvalue λ max B , the NB centralityũ in can be obtained through the much smaller matrix M [35] as shown in Appendix G.
Based on Eq. (44), the probability of newly infection k∈∂i αψ k→i (t) + βφ k→i (t) can be identified as the incoming vectors of the messages as αψ in i (t) + βφ in i (t), which will be increasingly more aligned withũ in i as time evolves. Thus, we can useũ in to predict the relative strengths of the outbreak, indicated by {1 − P i S (t)}. Figure 6 (a) demonstrates that, for β large enough, the evolution of correlation coefficient ρ betweenũ in and the profile of the outbreak on the WP2015 network generally increases with time. For small β = 0.002, the correlation remains low as the disease will die out. For a rather large β = 0.016, the correlation coefficient ρ increases rapidly in the initial stage of the development of the spreading, and then decreases to a lower level. This is because for a very large β, most nodes are likely to be infected eventually, irrespective of the spatial structure of the network. Such relation between ρ and β is also observed in other networks as shown in Fig. 6(b) .
Similarly, the IBMF approach dictates that the leading eigenvector of the adjacency matrix A (known as the eigenvector centrality [43] ) is a predictor of the outbreak profile [27] . Comparison between different centrality measures is briefly discussed in Appendix D, where the NB centrality generally provides a better prediction than the eigenvector centrality. As was mentioned before, the DMP approach only avoids the effect of mutual Figure 6 . Correlation coefficient ρ between the nonbacktracking centrality u in i and the outbreak profile, measured by the probability that each node i has been infected, i.e., 1 − P i S (t). The parameters are ν = 0.2, µ = 0.125. At time t = 0, there are 5 initial exposed nodes. (a) Correlation coefficient ρ of the spread on the WP2015 network as a function of time. (b) Correlation coefficient ρ as a function of β in the large time limit, where the spreading processes have saturated. Each data point is averaged over 5 instances with different sets of randomly selected initial exposed nodes.
infection due to one-step backtracking, and the approximation accuracy can deteriorate if counteracting the effect of one-step backtracking is insufficient. This is of particular concern when the NB centralityũ in displays the localization phenomenon, where the centrality values of a few nodes are much larger than the others [43, 44] . In Appendix D we demonstrated the degradation of the approximation power of the DMP equations and the NB centrality in random networks with a relatively large planted clique (i.e., a complete subgraph), which possesses the localization property [43] ; we also observed a poor approximation accuracy in one of the contact networks from the SocioPatterns data, recorded in a high school in 2013 (HS2013) [45] , where the NB centralitỹ u in is more localized on a few communities.
In addition to the social distancing rules that lower the transmission probabilities β and α, reducing social contacts between individuals is also an effective measure to slow down the spread of the disease. We examine its effects in the static contact networks considered here, by removing edges between nodes. Such measure will change the network structure and reduce the leading eigenvalue λ max B of the matrix B, which can enlarge the epidemic-free region in the parameter space as shown in Fig. 5(c) . In general, the leading eigenvalue λ max B can depend intricately on the network structure. In the special case of configuration model where the node degrees follow a given distribution P (d) and the nodes are wired randomly, it is found in [35, 43] that λ max B can be approximated as
which already indicates that even in the high degree limit, the epidemic thresholds obtained by R 0 in Eq. (43) do not generally coincide with those obtained in the DMP approach in Eq. (42) . Specifically, the second moment of the degree distribution is also relevant for epidemic thresholds in uncorrelated random networks [9] , which is not captured by Eq. (43) . A more refined approximation taking into account the relation of degrees of neighboring nodes is given in [44] . The accuracy of these approximations depends on the validity of the uncorrelated random network assumption and/or the presence of localization of NB centrality [44] (see also Appendice G 4). In Fig. 7(a) , Erdős-Rényi (ER) random graphs and networks generated via a stochastic block model (SBM) are considered. A network generated by SBM has 4 communities, where each community comprises 50 nodes; node i in community a is connected to c ab nodes of community b on average. Here, we consider c ab = 4, ∀a = b and c aa = 10, ∀a = 1, while different values of c 11 are considered. When c 11 = 10, all four communities are statistically equivalent, and the node degrees follow a Poisson distribution, similar to the ER random graphs. For such Poisson random graphs, λ max B ≈ d according to Eq. (46) , which is justified experimentally by the numerical results of Fig. 7(a) . On the other hand, λ max B deviates from d in SBM networks with c 11 > 0. The networks become sparser when edges are removed randomly, and λ max B decreases linearly with d . A quasi-linear decreasing trend is also observed in the random edge removal experiments in contact networks from the SocioPatterns collaboration as shown in Fig. 7(b) , as well as in scalefree networks where the node degrees follow the powerlaw distribution P (d) = d −γ / d max d=d min d −γ , as shown in Fig. 7(c) .
In Fig. 7(c) , the leading eigenvalues λ max B of scale-free networks can deviate significantly from d , especially for more heterogeneous networks with a small γ value. Therefore, predicting the cause of the spread through R 0 becomes very unreliable in scale-free networks, an effect which has been observed in network epidemiology studies [9] . Physically, there exist a small number of hubs (i.e., nodes with very large degrees) in scale-free networks, which can be viewed as super-spreaders that significantly facilitate the spread of the disease. In light of this, restricting contacts of these high-degree nodes preferentially can effectively reduce λ max B as shown in Fig. 7(d) , and consequently lower the epidemic thresholds.
In this work, we studied the SEIR model with presymptomatic transmissions, a feature of the COVID-19 disease, in both artificial and realistic contact networks through the dynamic message-passing method. The DMP approach provides a much better approximation compared to the IBMF approach, while being much less [46] , in a high school (HS2013) [45] and workplaces (WP2013, WP2015) [47, 48] . (c) Scale-free networks generated from the configuration model. The network size is N = 400 and the node degree is limited with d < d max = 100. (d) Scale-free network with γ = 2 (same as the one in panel (c)). In the edge removal process, a node i is firstly selected according to the probability pi = computational demanding than MC simulations. The linear stability analysis of the DMP equations gives rise to the epidemic thresholds and phase diagrams of the models, where their dependence on epidemiological parameters and the network structure are elucidated. A larger pre-symptomatic transmission probability value α leads to a lower critical point β c , which makes the strategy of blocking only symptomatic transmission less effective. We also show that different intervention strategies impact on the epidemic thresholds in different manners. The influence of network structure on the epidemic thresholds is represented by the leading eigenvalue λ max B of the non-backtracking matrix B, which encodes more subtle structural information in contrast to the average number of contacts d appeared in the basic reproduction number R 0 . Additionally, we demonstrated that the nonbacktracking centralityũ in related to the leading eigenvector of the matrix B can effectively predict the relative strength of the outbreak.
On the other hand, it is worthwhile mentioning some limitations of the DMP method. Firstly, as the DMP approach is based on the decorrelation assumption of the infection signals, it may become a less accurate approximation when the correlations between trajectories are non-negligible, which can happen when there is a single initial exposed node seeding the dynamics and/or there are many short loops in the network. Secondly, the approximation accuracy of the DMP approach also deteriorates when counteracting one-step backtracking reaction is insufficient to avoid the effect of mutual infection; this effect has been observed in networks where the non-backtracking centrality exhibits a localization phenomenon as shown in Appendix D. It is an interesting future direction to further characterize the condition of non-backtracking centrality localization, its impact on spreading processes and possible improvements [44] .
The theoretical frameworks were mostly applied to contact networks in some specific scenarios or those exhibiting particular characteristics, such as the presence of community structure or high-degree hubs. The applications in a wider scale (e.g., in a city) require considering additional network characteristics, such as the mixing pattern of different age groups [30] , the house-hold structure and so on [49] . Additional states, such as hospitalized and dead, can also be considered in order to model the pressure on public-health services and social cost. Since pre-symptomatic transmissions make it more difficult to contain the disease by dealing with the symptomatic cases only, an extension of particular interest is to examine the effectiveness and limitation of (manual or digital) contact tracing, mass testing and other strategies which can identify exposed individuals that have not shown symptoms. The DMP equations developed here will also benefit future works which aim at optimal deployment of resources (e.g., vaccines) to contain the spread of epidemics with pre-symptomatic transmissions [32, 50] . 
In this section, we derive the DMP equations from the principled dynamic belief propagation established in [16] . It is based on the message m i→j ( σ i | σ j ), which is the cavity probability of the dynamical trajectory σ i = [σ 0 i , ..., σ T i ] (where σ t i ∈ {S, E, I, R}) of node i in the cavity graph in which node j has been removed. Since the transition between states is irreversible (only permitted in the order S → E → I → R), the trajectory σ i can be parametrized by three transition times (τ i , ω i , ε i ) as σ i = |S 0 SSE τi EEI ωi IIIR εi RR T .
The dynamic belief propagation for the modified SEIR model takes the following form
where W SEIR is the transition kernel
The marginal of a trajectory of node i is computed as
The cavity probability of a trajectory has the following properties
where similar relations hold for m i (τ i , ω i , ε i |·).
The cavity probability of node i being in a certain state at time t is obtained by tracing over the cavity probability of trajectories,
where similar relations holds between P i σ (t) (σ ∈ {S, E, I, R}) and m i (τ i , ω i , ε i |·). Using the above definition, we compute the cavity probability P i→j S (t + 1) as
where we have assumed τ j = T in the cavity graph and traced over the τ i , ω i , ε i [16] ; the message θ k→i (t + 1) is defined as
which has the physical meaning of the cavity probability that node k (in either exposed or infected state) has not transmitted the infection signal to node i up to time t + 1. In order to eliminate the explicit dependence on the microscopic trajectories, we compute the iteration scheme of θ k→i (t + 1) as
where we have introduced the message φ k→i (t) (the cavity probability that k is in state I but has not transmitted the infection signal) and ψ k→i (t) (the cavity probability that k is in state E but has not transmitted the infection signal). The message φ k→i (t) is computed as
To proceed further, we first compute the update rule of ψ k→i (t) as
We observed that the message ψ k→i (t) can also be expressed in another form
Comparing Eq. (A17) and Eq. (A18) yields
Inserting Eq. (A19) into Eq. (A16) gives the update rule of φ k→i (t)
which closes the update rules of the messages θ i→j , φ i→j , ψ i→j and P i→j S . Collecting the incoming messages to node i gives rise to the marginal probability of node i being in state S
3. Probabilities of Being in State E, I, R
Finally, we consider the cavity probabilities P i→j E (t), P i→j I (t) and P i→j R (t). We first compute
and notice that
where we have made use of the property Eq. (A4). Inserting Eq. (A24) into Eq. (A22) gives rise to
Similarly, we compute
and make use of the property Eq. (A6) to compute
Inserting Eq. (A28) into Eq. (A26) gives rise to
Upon obtaining P i→j S (t+1), P i→j I (t+1), P i→j R (t+1), the cavity probability P i→j E (t+1) is given by the normalization condition
Applying the same analysis to the marginal of trajectory m i (τ i , ω i , ε i |·) leads to similar relations of the probability of node i being in state E, I, R as
which closes the DMP equations.
As for the COVID-19 epidemic, there is a relatively small number of patients who remain asymptomatic prior to recovery. Here, we discuss how to accommodate potential disease transmissions due to having asymptomatic individuals in our framework. Patients who do not exhibit symptoms at the time of exposure but developed symptoms at a later time are termed pre-symptomatic as described in the main text. There are many possible scenarios to model asymptomatic infections. One example is that each individual has a certain probability to become either asymptomatic or pre-symptomatic upon contracting the virus, in which case an additional asymptomatic state is required to model such stochastic transitions. The DMP equations can be derived accordingly, but will differ from those in this study.
Another perspective is based on the observation that whether an exposed individual will develop symptoms not seems to vary from person to person, based on age and pre-existing medical condition, e.g., children are more likely to have mild or no symptoms [23] . In light of this, one can assign each node i a label ℓ i (either probabilistically or according to additional existing information information), such that ℓ i = 1 if the individual is expected to develop symptoms when exposed and ℓ i = 0 otherwise. Nodes with different labels may have different epidemiological parameters. The disease transmission dynamics still obeys the transition rules of Sec. II, and the DMP equations in Sec. III B also apply. The only difference is that when a node i is in state I, one should interpret it as infected with symptoms if ℓ i = 1, or asymptomatic if ℓ i = 0. In this way, our theoretical framework can readily accommodate asymptomatic transmission, except that a weighted-version of non-backtracking is needed to accommodate different characteristics of symptomatic and asymptomatic individuals. The realistic contact networks are taken from data sets obtained from the SocioPatterns collaboration website [31], where the face-to-face contacts were recorded through wearable sensors over a certain period. Each data set contains a lists of active contacts between two individuals lasting for 20 seconds and the membership information of each individual (belonging to a class or department). To build the contact network, we first aggregate the contacts between any two individuals i and j, and consider the link between node i and node j as active if the cumulative contact duration between them in the recording period is not less than 60 seconds. A more precise treatment is to retain the information of contact duration in the transmission probabilities β ji , α ji . However, this results in a graph with a weighted non-backtracking (NB) matrix, which complicates the analysis. For simplicity, we only preserve the topological information of the resulting networks, keeping in mind that the contact duration information can also be incorporated in our framework. In this work, we consider contact networks in a primary school (PS2014) [46] , a high school (HS2013) [45] and workplaces (WP2013, WP2015) [47, 48] .
Artificially generated networks are also considered in this work, including Erdős-Rényi (ER) random graphs, random networks with community structures, scale-free networks and random networks with planted subgraph structures.
Random networks with community structures are generated through the stochastic block model. It is specified by a list of population of each block (n 1 , n 2 , ..., n r ) and a r × r symmetric probability matrix P ab , where r is the number of blocks. Suppose node i is assigned to block a and node j is assigned to block b, then node i and node j are connected with probability P ab . The average number of neighbors of node i (assigned to block a) belonging to block b is given by
The scale-free networks are generated through the configuration model. We first generate a degree sequence (d 1 , d 2 , ..., d N ) of size N where each element follows the power law distribution
where d min and d max are the minimum and maximum of the admissible degrees. In this study, we set d min = 5 and d max = 100, considering the maximal number of people that one contacts can not be arbitrarily large. After assigning a degree to each node, we randomly connect different nodes such that each node i has d i connections. The resulting graph may have a few self-loops and multiple edges between two nodes, which are simply removed to form a simple graph.
We also look at random networks with planted subgraph structures, primarily for the purpose of examining the effect of localization. Following [43] , we consider the ER random graph (with average degree d ) with a planted hub (of degree d h ), which is constructed by adding a hub to the existing ER network through creating connections from the hub to d h randomly selected nodes. When d h is large enough, i.e., d h > d ( d + 1), it is argued that the eigenvector centrality (i.e., the leading eigenvector of the adjacency matrix) is localized at the hub, while the non-backtracking Figure 8 . Evolution of the average number of individuals in state E or I, i.e., i (P i E (t) + P i I (t)) vs t. The underlying network is a random regular graph with N = 100, d = 10. The parameters are T = 100, ν = 0.2, µ = 0.125, β = 0.03, α = β/2. The initial exposed nodes are randomly selected. (a) One initial exposed node. (b) Two initial exposed nodes. (a) Five initial exposed nodes. In general, when there are more initial seeds, the approximation accuracy of the theory becomes better. centrality does not suffer from localization [43] . The intuition behind the localization of eigenvector centrality is that the hub and its neighbors are reinforcing each other, similar to the mutual infection effect of the IBMF approach to epidemic spreading. The NB centrality avoids this problem due to the hub by forbidding one-step backtracking.
On the other hand, it is noticed in [43] that a relatively large clique (i.e., a complete subgraph) can cause the NB centrality to localize as well. Therefore, we consider the ER random graph (with average degree d ) with a planted clique (of size N clique ), which is constructed by randomly selecting N clique nodes of the existing ER network to form a complete subgraph. For the original ER graph, the leading eigenvalue of the NB matrix B is given by λ max
In the presence of the clique, the leading eigenvalue satisfies λ max B > N clique − 2, which can easily exceed d for large N clique , so that the NB centrality is dictated by the clique [43] . There are also other subgraph structures that can cause the NB centrality to localize, such as dense subgraphs and overlapping hubs [44] . The intuition behind the localization of the NB centrality is that there is a subgraph sharing many neighbors, and avoiding one-step backtracking is insufficient to counteract the self-reinforcement among them. In this section, we give more examples of spreading processes experiments to support the findings in the main text. In Fig. 8 , it is shown that the approximation accuracy of the theories improves as the number of initial exposed nodes increases; we expect that it also depends on the locations of the initial seeds.
In Fig. 9 , we examine the theoretical results on ER random networks with different average degrees and/or different planted subgraph structures. We remark that there is no absolute fair comparison among different networks, as they have different dependencies on the epidemiological parameters (e.g., as in Fig. 3 of the main text) . Here, we fix the values of ν, µ and let α = β/2, and choose β such that approximately 70% of the population have contracted the disease in the final time T = 100. Fig. 9 demonstrates that the IBMF approach becomes a better approximation when the network is denser. This implies that the DMP approach is superior to the IBMF approach, especially in sparsely connected networks. When the average degree becomes higher, the trajectories obtained by the IBMF equations is approaching those by the DMP equations. In the limit of very dense random network, the mass-action approximation becomes a good approximation [9] . As mentioned before, a planted hub of relatively large degree can cause the leading eigenvector of the adjacency matrix A to localize [43] , which impairs the accuracy of the IBMF approach. We show in Fig. 9 (b) that such a planted hub does not have a noticeable effect on the approximation accuracy of the DMP approach. However, a planted clique which can cause the non-backtracking centrality to localize [43, 44] , does impair the approximation accuracy of the DMP approach as shown in Fig. 9 (c).
In Fig. 9 , we examine the theories on the networks extracted from contact data obtained in the SocioPatterns collaboration. The approximation accuracy of the HS2013 network is much poorer than the other networks, which may be attributed to the weakly localization of the non-backtracking centrality as shown in Fig. 13 . IBMF DMP MC (f)
The networks are ER random graphs of size N = 200 with different average degrees and/or different planted subgraph structures. The systems start with 5 initial exposed nodes. The parameters are T = 100, ν = 0.2, µ = 0.125, α = β/2, while β is selected such that in the final time T = 100, there are approximately 70% of the population that have contracted the disease. Comparing panels (a)(d)(e)(f), it is observed that the IBMF approach becomes a better approximation when the network is denser (i.e., d is larger). In panel (b), a planted hub of degree d h = 40 is created in the networks; the DMP approach is still a rather good approximation in this case. In panel (c), a planted clique (or complete subgraph) of size N clique = 20 is created in the networks; the approximation accuracy of the DMP approach is comparable poorer than other networks.
In Fig. 4 of the main text, the normalized outbreak size at large time r = i P i R (∞)/N was used to identify the critical point β c of the SEIR model. In Fig. 11(a) , we show the transient evolution of the normalized outbreak size, which is defined by the fraction of nodes that has contracted the disease as 1 − i P i S (t)/N . The normalized outbreak size increases only gently near criticality, while it grows rapidly above the critical point β c .
As the epidemic outbreaks are triggered by a few initial seeds, there is always a small probability that the disease will die out before spreading out further (e.g., the infected node transforming into state R much faster than average, or the infection signal not being transmitted to neighbors). This can also happen when β > β c , indicating that there are some instances where the outbreak size is small although the system is in the global-epidemic phase, as shown in Fig. 11(b) . This phenomenon has been observed in the SIR model [28, 39] . While the theoretical frameworks only concern the average behaviors, they do not capture such variability of trajectories due to stochastic fluctuations. A detailed theoretical investigation into these aspects will be an interesting topic for future studies.
As mentioned in the main text and in Sec. D 1, the degradation in accuracy of the approximation of the theory (IBMF or DMP method) is correlated with the localization phenomenon of the corresponding centrality measure, where the centrality values of a few nodes are much larger than the others. We have showcased this in Fig. 9 by considering some planted subgraph structures in a random network. We further examine the prediction of the centrality measures on the outbreak profiles in these planted random networks. In Fig. 12(a) , it is shown that in the ER graph, both the eigenvector centrality and the NB centrality (coming from the linear approximation of the IBMF and DMP approaches), as well as the degree, are good predictors of the outbreak profile, despite the poor approximation of the full nonlinear IBMF approach in Fig. 9(a) . On the other hand, a planted hub causing the eigenvecter centrality to localized, degrades its prediction accuracy as shown in Fig. 12(b) , while the NB centrality appears to be a much better predictor. In the presence of a clique, the NB centrality also predicts the outbreak poorly as shown in Fig. 12 (c) due to the localization phenomenon. Figure 10 .
The networks are extracted from contact data obtained in the SocioPatterns collaboration. The systems start with 5 initial exposed nodes. The parameters are T = 100, ν = 0.2, µ = 0.125, α = β/2, while β is selected such that in the final time T = 100, there are approximately 70% of the population that have contracted the disease. The approximation accuracy of the HS2013 network is much poorer than the other networks. In all three cases, the degree appears to be a good predictor of the outbreak profile; in random networks, having more neighbors usually implies a higher chance to contract the disease. However, this does not hold in general, e.g., consider a hub connected to many dangling nodes but linked to the bulk of the network through only a few edges, such a hub node is not likely to be at high risk as indicated by its degree.
It has been shown in Fig. 10 that the approximation accuracy of the HS2013 network is much poorer than the other networks from the SocioPatterns data sets. In Fig. 13 , we show the NB centralities on the corresponding networks. Comparing to other networks, the NB centralityũ in of HS2013 concentrates on two communities, which may cause the approximation of the DMP equations to be less accurate. eigenvector can be deduced if the matrix X is irreducible [37] .
In the main text, we have shown that the eigenvalue of the Jacobian J and the eigenvalue of the non-backtracking matrix has the relation
which can be solved for λ J , leading to
From this relation, one can figure out the principal eigenvalue of J . Noticing that λ max
= ρ(B) ≥ 0 due to the Perron-Frobenius theorem, we can focus on λ + J and λ B ∈ R since the eigenmode with a fastest growth rate will not be realized by negative or complex eigenvalue of J . We also assume that ρ(J ) and ρ(B) are non-zero. To simplify the notation, denote a = (
where we have made use of the fact that 2bβ α > 0 and (a−c+αλB )
The epidemic threshold is obtained by solving λ max J (β, α, ν, µ, λ max B ) = 1.
Similar to the DMP approach, we can also derive the epidemic thresholds through the individual-based mean-field (IBMF) approach. The initial disease-free state is perturbed infinitesimally as P i S (0) = 1 − ǫ i , in which case the probabilities P i E (t) and P i I (t) are also of order ǫ i in the initial stage. We expand Eq. (8) and neglect terms of higher order of ǫ i , leading to
Equations (E5) and (9) constitute a linear dynamical system of the probabilities {P i E (t), P i I (t)}. They can be written in the matrix form as
where the 2N × 2N Jacobian matrix J MF is defined as
where I is the identity matrix and A is the adjacency matrix of the graph. The spectral radius ρ(J MF ) determines the growth rate of the fastest mode of Eq. (E6). Due to the Perron-Frobenius theorem, ρ(J MF ) equals the leading eigenvalue λ max J MF of J MF , which is related to the leading eigenvalue λ max For completeness, we provide a derivation using the matrix notation in this section. To this end, we define the N × 2|E| source matrix S and target matrix T S ie = 1 if node i is the source of directed edge e, i.e., ∃j, e = i → j, 0 otherwise, (G2)
T ie = 1 if node i is the target of directed edge e, i.e., ∃j, e = j → i, 0 otherwise,
and a 2|E| × 2|E| auxiliary matrix
where I 2|E| is the 2|E|-dimensional identity matrix. We further index the set of directed edges according to the order (e 1 , ...e |E| , e |E|+1 , ...e 2|E| ) such that for n ≤ |E|, e n = i → j, one has e n+|E| = j → i. Under this notation, the non-backtracking matrix B can be written as
It can also be easily verified that
The key element of the derivation in Ref. [35] is that for a given 2|E|-dimensional u, one defines the corresponding N -dimensional incoming and outgoing vectors 
where we have made used of the relations in Equations (G6) and (G7). The above expressions can be written as
where we identify the matrix M defined in Eq. (G1). Now suppose u is an eigenvector of B with eigenvalue λ B , then Eq. (G12) leads to
which implies that Su T u or u out u in is the eigenvector of M with eigenvalue λ B . Therefore, we can work with the much smaller matrix M for computing the spectrum of the non-backtracking matrix B.
Physically, this reduction comes from compressing the edge-based data {u i→j } to node-based data {u out i , u in i }. In the context of spreading processes, they correspond to edge-based messages and node-based marginal probabilities, e.g., the cavity probability P i→j S in the linearized dynamics satisfies P i→j S (t + 1) ≈ P i→j S (t) − k∈∂i\j αψ k→i (t) + βφ k→i (t) ,
while the marginal probability satisfies
where the probability of newly infection k∈∂i αψ k→i (t) + βφ k→i (t) corresponds to the incoming vectors of the messages as αψ in (t) + βφ in (t).
or explicitly
The above equations can be reduced to a nonlinear eigenvalue problem
and implies a relation between the outgoing-component u out and the incoming-component u in of the eigenmode as
In the main text, it is argued that the leading eigenvector of the matrix B (denoted asũ), and the corresponding incoming vectorũ in = Tũ (i.e., the non-backtracking centrality), are useful in predicting the outcome of the epidemics. From the analysis in this section, both the leading eigenvalue λ max B and the non-backtracking centralityũ in can be obtained through the matrix M , which significantly reduces the computational complexity. The Perron-Frobenius theorem guarantees that the leading eigenvectorũ of the matrix B can be chosen to be non-negative, so does the non-backtracking centralityũ in .
The leading eigenvalue of the matrix M can be computed exactly for random regular graphs. We first notice that the smallest eigenvalue of the Laplacian matrix L = D − A of a graph is λ min L = 0 [52] . For regular graphs with degree d, we have D = dI, and therefore L commutes with A. It suggests that the largest eigenvalue of A is
Equation (G18) can be expressed as
which implies that λ M + d−1 λM is an eigenvalue of the adjacency matrix A with eigenvalue u in , denoted as λ A . The eigenvalue λ M of M is related to λ A as
Therefore, the leading eigenvalue λ max M of M , which is equal to the leading eigenvalue λ max B , is identified as can be derived, which can be useful for estimating epidemics properties in large networks. In the annealed approximation where only the information of the degree distribution P (d) is retained, λ max B is approximated as [35, 43] λ max,an
In a more refined approximation assuming uncorrelated networks but taking into neighborhood information, λ max B is approximated as [44] λ max,un
.
To examine the validity of the uncorrelated random network assumption, we consider the degree correlation coefficient [53] 
kl(e kl − q k q l ),
where e kl is the element of the degree correlation matrix (i.e., the probability of finding an edge connecting two nodes of degree k and degree l), q k = kP (k) d (with P (k) being the degree distribution) is the probability of finding an edge whose one end has degree k, and σ 2 = d max k=1 k 2 q k − d max k=1 k 2 q k 2 is the variance of the measure q. The neutrality condition r d = 0 needs to be (at least roughly) satisfied for a network to be uncorrelated. In Table. I, we examine the approximations offered by Eqs (G24) and (G25) for the networks considered in this work. In general, λ max,un B provides a better approximation than λ max,an B , both of which predict λ max B quite well when r d is low. The two cases with a relatively poor approximation (especially for λ max,an B ) are the SBM network with c 11 = 30 and the ER network with a clique. These two networks exhibit high values of the degree correlation coefficient r d , violating the assumption of uncorrelated random network. On the other hand, the existence of dense subgraph structures can cause localization of the NB centrality, which also makes the approximation inaccurate. In light of this, it has been proposed in Ref. [44] to identify some characteristic subgraph structures for an improved approximation of λ max B . 