key: cord-024746-ijlnefz3 authors: Koher, Andreas; Lentz, Hartmut H. K.; Gleeson, James P.; Hövel, Philipp title: Contact-Based Model for Epidemic Spreading on Temporal Networks date: 2019-08-02 journal: nan DOI: 10.1103/physrevx.9.031017 sha: doc_id: 24746 cord_uid: ijlnefz3 We present a contact-based model to study the spreading of epidemics by means of extending the dynamic message-passing approach to temporal networks. The shift in perspective from node- to edge-centric quantities enables accurate modeling of Markovian susceptible-infected-recovered outbreaks on time-varying trees, i.e., temporal networks with a loop-free underlying topology. On arbitrary graphs, the proposed contact-based model incorporates potential structural and temporal heterogeneities of the contact network and improves analytic estimations with respect to the individual-based (node-centric) approach at a low computational and conceptual cost. Within this new framework, we derive an analytical expression for the epidemic threshold on temporal networks and demonstrate the feasibility of this method on empirical data. Accurate models of disease progression are valuable tools for public health institutions as they enable detection of outbreak origins [1] [2] [3] [4] , assessment of epidemic risk and vulnerability [5] [6] [7] [8] , and containment of the spreading at an early stage [5, 9] . Mitigation strategies can thus be evaluated and employed without the need to run a large number of Monte Carlo (MC) realizations. A fundamental challenge to mathematical epidemiologists is the accurate determination of the critical parameters that separate local and global epidemic outbreaks [10] [11] [12] [13] [14] [15] [16] . To this end, the early Kermack-McKendrick model [17] divides a population according to the disease status into compartments of susceptible, infected, and recovered individuals with mass-action equations to determine the transitions between them. Since then, a wide range of improvements has been proposed, including the impact of stochasticity [18] [19] [20] , non-Markovian dynamics [21] [22] [23] [24] [25] , and, notably, heterogeneity in the contact structure [11, [26] [27] [28] [29] [30] . In recent years, the availability of mobility and contact data with a high temporal resolution, so-called temporal networks, offers another opportunity to improve analytical predictions [31] [32] [33] [34] [35] [36] . The timing of links between nodes matters, in particular, when the network evolves on a similar timescale as the spreading dynamics, which led to an increasing interest in the interplay between disease and network dynamics [37] [38] [39] [40] [41] [42] . One approach to model the states of individual nodes in a network takes the corresponding probabilities directly as variables in a set of coupled dynamic equations [7, 13, 20, [43] [44] [45] [46] [47] . We refer to this approach as the individualbased (IB) model, though it is sometimes also called the N-intertwined model [20] or quenched mean field [48, 49] . However intuitive, the analytic predictability suffers from the simplifying assumption that epidemic states of adjacent nodes are independent. Recently, a change from a node-centric to an edgecentric perspective has been discussed within different frameworks in order to overcome the inherent limitation of the IB model. These approaches include branching processes [50] , message passing [23, 51] , belief propagation [4] , and the edge-based compartmental model [39] . So far, however, edge-centric models are mostly limited to static topologies. It thus remains an open challenge to simultaneously account for topological and temporal properties of the underlying contact data and hence improve current predictions of the epidemic threshold [7, [52] [53] [54] [55] [56] . In this paper, we generalize the dynamic message-passing approach for discrete-time Markovian susceptible-infectedrecovered (SIR) spreading [51] to time-evolving networks and derive the epidemic threshold within this new framework. The proposed model takes an edge-centric perspective because the relevant dynamic equations are based on the set of edges. Furthermore, the framework integrates the complete temporal and topological information of the underlying network into the epidemic model. We refer to our approach as contact-based (CB) model and compare numerical predictions with the widely used IB model that takes a node-centric perspective. Within the CB framework, we then derive a new analytic expression of the epidemic threshold for temporal networks and show that the edge-centric approach improves existing results [7, 20, 43, 44, 52, 53] at a low conceptual and numerical cost. The CB and IB models have been implemented in PYTHON with the source code available on Github [57] . The remainder of this paper is structured as follows: First, we summarize the conceptual framework in Sec. II and formulate in Sec. III the dynamic equations of the IB and CB models. Then, we derive the epidemic threshold for temporal networks within the CB framework in Sec. IV. We compare the edge-and node-centric approaches against MC simulations in Sec. V and close with a discussion in Sec. VI. Appendix A includes an extension to weighted contacts and heterogeneous epidemiological parameters. A network analysis of the German cattle-trade data is given in Appendix B. Further results and applications of the CB model are summarized in Appendix C. We consider a temporal network G ¼ ½Gð0Þ; Gð1Þ; …; GðT − 1Þ with N nodes and T snapshots sampled at a constant rate. Although both modeling frameworks can, in principle, account for contact weights that indicate the strength of a connection, we focus on unweighted networks for simplicity and refer to Appendix A for an extension of the model. Emphasizing the important difference between temporal and static elements, we refer to contacts as time-stamped links ðt; k; lÞ ∈ C ⊂ T × N × N , thereby denoting with N , T , and C the set of nodes, time stamps, and contacts, respectively. We further assume that every contact is of constant duration and equal to the sampling time of the temporal network. By edges, we refer to the corresponding static elements in the time-aggregated network. In other words, an edge ðk; lÞ ∈ E ⊂ N × N exists if and only if at least one (temporal) contact is recorded between k and l. Here, we denote with E the set of edges. Moreover, we assume directed edges throughout the paper and represent an undirected contact as two reciprocal contacts. Following the convention in Ref. [23] , we denote with k → l a directed edge from k to l, and we indicate edge-based quantities in a similar fashion. As the stochastic process, we assume a discrete-time SIR model, where a node l ∈ N represents an individual that is either susceptible, infected, or recovered at a given time t with a corresponding probability S l ðtÞ, I l ðtÞ, and R l ðtÞ, respectively. A susceptible node that is in contact with an infected neighbor contracts the disease with a constant and uniform (per time step) probability β. Furthermore, we treat the transmission events from multiple infected neighbors as independent, and similarly, we interpret potential (integer) edge weights as independent infection attempts (see Appendix A). We do not account for secondary infections within one time step; i.e., only direct neighbors can be affected. Once infected, the individual recovers with a uniform and constant probability μ independently of the infection process and henceforth acquires a permanent immunity. Concerning the contact data, we focus our numerical analysis first on a face-to-face interaction network between 100 conference participants [58] . This so-called proximity graph has a resolution of 20 s, and the observation time is limited to the first 24 h. If necessary, we extend the data set with a periodic boundary condition in time. The timeresolved contacts enable the study of spreading of airborne diseases as well as the propagation of ideas and rumors. As an illustrative example, we present in Fig. 1 the timedependent probability that a selected node in the proximity graph is either susceptible (yellow), infected (red), or recovered (green). The results are derived from 10 4 MC simulations FIG. 1 . Illustrative examples of a simulated epidemic outbreak from a single initially infected node. Colors give the probability that another arbitrarily selected node is in the susceptible (yellow), infected (red), or recovered (green) state, respectively. Simulation parameters: μ ¼ 2.85 × 10 −4 , β ¼ 100μ, 10 4 MC realizations. with the same initially infected node. The trajectories reflect the bursty activity of the underlying temporal network [58] within the first 12 h and the subsequent inactive nighttime. As a second source of data with direct relevance to public health, we consider an excerpt of the national German livestock database HI-Tier [59] . This temporal network comprises the movement of cattle between farms in Germany for the year 2010 with daily resolution. Within the observation window of 365 days, more than 3 million transactions have been recorded between over 180 000 farms and traders, respectively. For more details on the graph, see Appendix B. Cattle trade is considered an important transmission route for livestock-related diseases such as foot-and-mouth disease (FMD), which broke out in the United Kingdom in 2001 with an estimated cost of 8 billion pounds sterling [60] . Therefore, the analysis of the corresponding spatiotemporal graphs is highly relevant to public health institutions. In this section, we present the mathematical framework to model the stochastic SIR process as outlined in the Introduction and Sec. II. Our main focus is the CB model, but in order to facilitate a direct comparison between the node and the edge-based approach, we begin with a short overview of the IB model. In the IB model, the marginal probabilities S l ðtÞ, I l ðtÞ, and R l ðtÞ for all l ∈ N directly enter a set of 3 × N coupled dynamic equations. The probability for l to contract the infection from k upon a temporal contact is given by βI k ðtÞ. For convenience, we introduce an indicator function with a k→l ðtÞ ¼ 1 if a (directed) contact from k to l exists at time t and a k→l ðtÞ ¼ 0 otherwise. Then, the probability for node l to receive no infection at time t from any of its neighbors factorizes by assumption to Q k ½1 − βa k→l ðtÞI k ðtÞ and k ∈ N . With this result, the marginal probability S l ðt þ 1Þ can be expressed by the probability S l ðtÞ to be susceptible in the previous time step t and not contract the infection within the interval ½t; t þ 1Þ. In the IB model, the joint probability factorizes, and we obtain Here, the crucial simplification is to treat the epidemic states of l and its neighbors as mutually independent, which is sometimes referred to as neglecting dynamic correlations [61] . The marginal probability I l ðt þ 1Þ follows from two independent contributions: (i) The outflux μI l ðtÞ indicates the transition from the infected to the recovered state. (ii) The influx ΔS l ðtÞ ¼ S l ðtÞ − S l ðt þ 1Þ reflects the probability that node l is newly infected at time t þ 1. Combining both contributions leads to The set of 2 × N coupled dynamic equations (1) and (2) thus constitutes the IB model for temporal networks. The remaining marginal probability R l ðtÞ to find node l in the recovered state follows from the conservation condition S l ðtÞ þ I l ðtÞ þ R l ðtÞ ¼ 1 for all l ∈ N . Finally, we assign a probability z l ¼ S l ð0Þ that node l is initially susceptible, as well as I l ð0Þ ¼ 1 − z l and R l ð0Þ ¼ 0 throughout the paper. Though intuitive and, in many cases, sufficient from a modeling perspective, the limits of the IB model are difficult to estimate due to the ad hoc factorization of the joint probability in Eq. (1). Even for the simplest network with two nodes connected by an undirected static edge, the IB approach can deviate significantly from the expected outcome as illustrated in Ref. [62] . In their example, recovery is neglected for simplicity, and only the first node is infected initially with some probability 0 < z 1 ≤ 1. Counterintuitively, the probabilities to find each node in the infected state converge to I 1 ð∞Þ ¼ I 2 ð∞Þ ¼ 1 according to the IB model, independent of the initial condition z 1 . This convergence occurs because integrating Eqs. (1) and (2) admits a probability flux from the outbreak location to the adjacent node and back to its origin again. This mutual reinfection, coined the echo chamber effect in Ref. [62] , appears because we neglect the fact that the probability I 2 to find the second node in the infected state is conditioned on the state of the first node, and thus the factorization in Eq. (1) is not justifiable. In an arbitrary network, an initially infected node leads to a cascade of secondary infections within which all marginal probabilities are highly correlated. An accurate model excludes these previously infected nodes from those that can potentially contract the infection in the future. In the next section, we discuss how a shift from a node-centric to an edge-centric view can take into account some such dependencies. We begin with a slightly different approach to the marginal probability S l ðtÞ. First, we note that l is susceptible at time t, if it was susceptible initially [with probability S l ð0Þ ¼ z l ] and has not contracted the infection from any of its neighbors up to time t. We assign the probability Φ l ðtÞ to the latter statement. Thus, without introducing any approximation at this stage, we can write In order to determine Φ l ðtÞ, we make the assumption that the underlying time-aggregated graph is a tree (ignoring directionality). Then, different branches originating in node l are independent as long as l remains susceptible, and thus Φ l ðtÞ factorizes. However, if node l contracts a disease from a neighbor k with some probability and passes it on to another node k 0 , then the corresponding probabilities I k and I k 0 are clearly correlated. A simple solution that allows different branches to nonetheless be treated as independent is to prevent a probability flow through the root node in the first place. From a graph-theoretic perspective, this solution corresponds to the (virtual) removal of all out-directed contacts from the root node. This approach does not modify the dynamics of the node under consideration because it can still contract the disease, and once infected, the recovery process is independent of the topology. However, the idea considerably reduces the amount of bookkeeping that would otherwise be necessary if we accounted for the correlations directly. The singular node l is said to be a cavity node or in the cavity state [23, 51] , a concept closely related to the test-node assumption [39] and the idea of cut vertices [63] . With this concept, we can factorize Ψ l ðtÞ and thus obtain Here, we introduce the probability θ k→l ðtÞ that no disease has been transmitted from node k to the cavity node l up to time t. The change in perspective towards an edge-centric analysis introduces new auxiliary dynamic quantities such as θ k→l ðtÞ. These quantities are defined on the set of edges E of the time-aggregated network, and thus the number of dynamic variables scales with L, the number of edges. In order to obtain a system of dynamic equations, we focus on our first edge-centric variable θ k→l . Initially, no disease was transmitted such that θ k→l ð0Þ ¼ 1 for all edges ðk; lÞ ∈ E. Henceforth, the dynamic quantity reduces only (i) upon a temporal contact indicated by a k→l ðtÞ and (ii) if the adjacent node k is infected without having transmitted the disease earlier to the cavity node l-we denote the corresponding probability by I k→l ðtÞ. Hence, the out-flow of probability is given by βa k→l ðtÞI k→l ðtÞ, leading to our first dynamic equation Next, the probability I k→l ðtÞ evolves according to three contributions. (i) It decreases with the recovery probability μ and (ii) with the probability β to infect its target node upon a temporal contact. These processes are independent and may contribute simultaneously with the joint probability βμ. (iii) I k→l ðtÞ increases with the probability ΔS k→l ðtÞ ¼ S k→l ðtÞ − S k→l ðt þ 1Þ that k is newly infected by at least one of its incident neighbors excluding the cavity node l. In sum and with the initial condition I k→l ð0Þ ¼ 1 − z k , these contributions lead to Finally, we consider the probability S k→l ðtÞ that node k, adjacent to the cavity node l, is susceptible. Since k is not affected by the state of l, it stays susceptible if it does not contract the disease from any of its remaining, incident neighbors j ∈ N k nl. It has been shown in Ref. [64] that the corresponding probability Φ k→l ðtÞ ¼ Q j∈N k nl θ j→k ðtÞ factorizes, and thus, similar to Eq. (3), we find S k→l ðtÞ ¼ z l Φ k→l ðtÞ or, equivalently, The disease progression in the CB framework is fully characterized by Eqs. (5) and (6), a set of 2L coupled equations. Equation (7) is introduced here for convenience only and can be substituted into Eq. (6). Next, we return to the node-centric quantities. To this end, we note that S l ðtÞ has already been determined in Eq. (4). The remaining marginals I l and R l are equivalent to the IB model and given by the conservation condition, as well as the transition to the recovered state in Eqs. (8) and (9), respectively: The CB model is exact for temporal networks, where the undirected, time-aggregated graph has a tree structure and is therefore loop-free. Most realistic networks, however, contain a large number of loops such as triangles in social graphs, where two friends are likely to have many more friends in common. Here, the CB model nevertheless appears to be "unreasonably effective" (cf. Ref. [65] ) and improves predictions significantly with respect to the IB approach as we will see in Sec. V. For further extensions to the model that include heterogeneous infection and recovery probabilities, as well as weighted contacts, see Appendix A. The parameters that mark the epidemic threshold can be derived by examining small perturbations around the disease-free state. If such perturbations die out, then any outbreak remains local, but if the perturbation grows, then a global epidemic may occur. We consider a linearization of the dynamic equations (5)-(7), which will give rise to a criticality condition, determining the epidemic threshold. We begin with the ansatz θ k→l ðtÞ ¼ 1 − δ k→l ðtÞ and z l ¼ 1 − ϵ l , where δ k→l ðtÞ; ϵ l ≪ 1 are small perturbations around the disease-free state for all nodes l and edges ðk; lÞ. Thus, Eq. (5) becomes In Eq. (7), we keep the linear terms of the Taylor expansion, which transforms the product into a corresponding sum: In Eq. (11b), we substituted the dynamic equation (10) and identified S k→l ðtÞ in the next step. From the resulting Eq. (11c), we can read the linearized form of ΔS k→l , which allows us to decouple the dynamic equations for I k→l : Next, we rewrite the remaining set of L dynamic equations in a compact, matrix-based formulation and therefore introduce the vectors IðtÞ and aðtÞ with elements I k→l ðtÞ and a k→l ðtÞ, respectively. To this end, we also express the linear operation P j∈N k nl a j→k ðtÞ in Eq. (12), which acts on the elements I k→l ðtÞ of the state vector, through the temporal unweighted nonbacktracking matrix BðtÞ: In other words, B k→l;j→k 0 ðtÞ ¼ 1 if the contact ðt; j; k 0 Þ at time t is incident on the edge ðk; lÞ (implying k 0 ¼ k), and additionally j ≠ l. Otherwise, we have B k→l;j→k 0 ðtÞ ¼ 0. It is only the nonbacktracking property j ≠ l that sets B apart from the adjacency matrix of the ordinary line graph. For temporal networks, a subtle distinction has to be made between the first and the second index of the L × L dimensional matrix B: The first corresponds to an outdirected (static) edge ðk; lÞ ∈ E of the underlying aggregated network and can be interpreted as a potential contact in the future. The second, however, is an incident (temporal) contact ðt; j; k 0 Þ ∈ C from node j to k 0 at time t. We also introduce the diagonal matrix diag(1 − βaðtÞ), with diagonal elements given by the vector 1 − βaðtÞ. Here, we denote by 1 the vector of all ones. With these definitions, we rewrite Eq. (12) as The explicit solution to the state vector IðTÞ at final observation time T is formally given by IðTÞ ¼ Pðβ; μÞIð0Þ, where the so-called infection propagator P [54] is introduced for notational convenience: In order to evaluate the asymptotic behavior, we assume a periodic boundary condition in time, i.e., BðtÞ ¼ Bðt þ TÞ. This allows us to assess the vulnerability of the temporal network through the spectral radius of the propagator P. In particular, we find that a SIR-type outbreak is asymptotically stable under small perturbations, i.e., remains confined to a small set of nodes, as long as the spectral radius satisfies ρ½Pðβ; μÞ < 1. Thus, the phase transition is given by the criticality condition Note that for irreducible and non-negative matrices, the largest eigenvalue is simple and positive according to the Perron-Frobenius theorem [66] . Assuming 0 ≤ β, μ < 1, a sufficient condition for temporal networks is to restrict contacts to the giant strongly connected component (GSCC) of the underlying time-aggregated graph. In Sec. V B, we fix the recovery probability μ and determine the critical infection probability β crit as the root of fðβÞ ¼ 1 − ρ½Pðβ; μÞ for different empirical networks. We conclude this section with a discussion on the static network limit. In the so-called quenched regime, the disease evolves on a much faster timescale than the dynamic topology and thus operates on an effectively static network with BðtÞ ≡ Bð0Þ ≡ B and aðtÞ ≡ 1 for all times t. As in the temporal analysis, we restrict the network to the GSCC so that the Perron-Frobenius theorem [66] applies. In this limit, the dynamic equations (5)-(7) reduce to the dynamic message-passing formulation in Ref. [51] . Moreover, Eq. (15) now becomes a product Q T−1 t¼0 P fast ðβ; μÞ ¼ ½P fast ðβ; μÞ T of T identical, single time-step propagators where 1 ¼ diagð1Þ denotes the identity matrix. The spectral radius in Eq. (16) factorizes to ρ½P fast ðβ; μÞ T ¼ ρ½P fast ðβ; μÞ T , and it follows that the criticality condition Eq. (16) reduces to ρ½P fast ðβ; μÞ ¼ 1. Furthermore, we find from basic linear algebra that ρ½P fast ðβ; μÞ ¼ ð1 − μÞð1 − βÞ þ βρðBÞ, and hence we obtain the corresponding static threshold condition β β þ μ − βμ The criticality condition in Eq. (18) deviates from the continuous-time result in Refs. [16, 23] . In the derivation presented here, the term βμ in Eq. (18) accounts for the simultaneous events when a node infects a neighbor and recovers within the same time step. In contrast to the quenched regime, one can also consider the so-called annealed limit. Then, parameters β and μ are sufficiently small such that no more than one infection or recovery event can take place within the observation time. Therefore, we expand the infection propagator to the first order in β and μ and obtain Here,ā ¼ 1=T P t aðtÞ andB ¼ 1=T P t BðtÞ denote the corresponding time-averaged quantities. It is insightful to evaluate simple bounds for the set of parameters ðβ; μÞ crit;slow that satisfy the threshold condition ρðP slow Þ ¼ 1 in the annealed limit. With 1=T ≤ā ≤ 1 for all elements inā, we thus find β Assuming the upper bound in Eq. (20) overestimates the outbreak risk and can be considered a conservative choice from an epidemiological perspective. This limit is realized for a temporal network where every edge appears exactly once within the observation time, henceā ¼ 1=T. The lower bound in Eq. (20) is exact in the case of a static network (thusā ¼ 1) and corresponds to the continuoustime result in Ref. [16] . However, this limit underestimates the outbreak risk, and therefore we conclude with a note of caution when applying results from static network theory directly to time-varying topologies. A big advantage common to both the node-centric IB and edge-centric CB modeling framework is a significant reduction in computational complexity compared to MC simulations. The CB model requires iteration through all edges at every time step, and thus the time complexity scales with OðLTÞ. The IB formulation and a single MC realization require OðCTÞ, whereC denotes the average number of active contacts, which can be significantly smaller than L. Stochastic MC simulations, on the other hand, require a large number of realizations in order to provide reliable statistics. The computational disadvantage of MC simulations becomes even more apparent when we consider a complex quantity such as the epidemic threshold, which requires multiple ensemble averages for different sets of epidemic parameters in order to fit the critical infection probability (see Sec. V B). Equally important, however, is the accuracy of our analytic approach. Therefore, in this section, we compare estimations from the IB and CB mean-field model with MC simulations using empirical data as introduced in Sec. II. We begin with an analysis on the level of individual nodes. In Fig. 2 , we show the cumulative infection probability for a small number of example nodes from the conference data set given the same outbreak location. The selection is intended to present qualitatively different trajectories, also demonstrating that deviations between the two models vary considerably. The MC result (blue curve) in Fig. 2 (a) corresponds to the introductory example in Fig. 1 . Here, a comparison with the analytic estimation shows that the CB approach leads to a substantial improvement to the IB model. Also in Figs. 2(b)-2(d), the trajectories are erratic, as they reflect the sudden changes in the underlying topology, highly individual and yet well approximated by the CB model. For all nodes in the network, we find that the CB model gives a closer upper bound to MC simulations because, unlike the IB framework, it accounts for dynamic correlations between nearest-neighbor states. Dynamic mean-field models such as the IB and CB framework provide realistic expectation values only if stochastic fluctuations are negligible. In order to illustrate the limitations, we study epidemic outbreaks for three different initially infected nodes in Figs. 3(a)-3(c), respectively. The left column gives the time-resolved distribution of the outbreak size, and the right column presents the final distribution at the end of the three-day observation period. For the ensemble average (blue line), we consider only realizations with more than 20 infected nodes overall. This threshold separates outbreaks that die out early due to stochastic fluctuations and thus permits a direct comparison with estimations from the IB and CB frameworks in green and red, respectively. We choose the outbreak locations such that the degree of stochasticity increases from top to bottom. In Fig. 3(a1) , we find a narrow distribution around the ensemble average, which is well approximated by the mean-field models. Minor outbreaks due to early extinctions are well separated in Fig. 3(a2) from large epidemics. In Fig. 3(b) , the initially infected node leads to realizations with considerably stronger fluctuations, and in Fig. 3(c) , it is barely possible to separate early extinctions at all. Additionally, we observe a second source of stochastic variation, namely, the time at which a disease escalates and hence evolves into a global epidemic. As a consequence, early outbreak sizes may be overestimated significantly before the analytic trajectory approaches the expectation value again [see Fig. 3(c1) ]. Remarkably, the performance of both mean-field models varies significantly with the outbreak location, even for the basic reproduction number R 0 well above the epidemic threshold. At the late phase of an outbreak, however, the mean-field models provide good approximations, and consistently with Fig. 2 , we find that the CB model outperforms the IB approach. In Appendix C, we demonstrate how a sufficiently large number of initially infected individuals significantly improves the predictability. Another source of stochasticity is the choice of disease parameters β and μ, respectively. We focus on the final outbreak size, averaged over all outbreak locations. The distribution as a function of the infection probability β (see Fig. 4 ) shows a percolation-like transition from localized spreading to epidemics that affect a considerable fraction of the network. We apply the same threshold as in Fig. 3 for a direct comparison between the averaged outbreak size and the mean-field models for β > 0.02. Here, we find that the difference between the expected size and the CB estimation is close to negligible, whereas the IB model consistently overestimates the expected value. A comparison at low values of the infection probability β becomes unreliable as stochasticity impedes a reasonable distinction between minor and global outbreaks. In order to illustrate the effect, we present in Fig. 5 the outbreak-size distribution for different values of β as marked by the arrows in Fig. 4 . This representation highlights the transition from the subcritical to the supercritical parameter domain: The unimodal distribution in Fig. 5(a) characterizes localized outbreaks, whereas the bimodal distribution in Fig. 5 (d) clearly separates early extinctions and global epidemics. Next, we focus on the critical infection probability that marks the transition. In Fig. 6(a) , we present the region of small β from Fig. 4 in order to focus on the transition from localized outbreaks to the sudden emergence of global epidemics. We determine the critical infection probability β crit (vertical blue line) from the maximum of the relative standard deviation [54] , also known as the coefficient of variation [see blue line in Fig. 6 Here, we denote with hσi and hσ 2 i the first and second moments of the outbreak-size distribution. The coefficient of variation captures the intuition that fluctuations dominate the outbreak-size distribution close to the transition. Indeed, c v diverges at the critical point for infinitely large networks, indicating a second-order phase transition [67] . Analytically, we determine β crit from the spectral criterion in Eq. (16) for the CB model and similarly within the IB framework [54, 68] . The comparison in Fig. 6 shows that the IB and CB models, marked by a red dashed and green dotted line, respectively, underestimate the critical infection probability from MC simulations (blue line) and thus overestimate the outbreak risk. Consistent with our previous results, we can state that a shift from a node-to an edge-centric framework improves the analytic prediction. In Appendix C, we present similar results for different values of the recovery probability μ. Next, we continue with a realistic application of the epidemic threshold to the German cattle-trade network. We now consider a completely different data set, where the system size is large and contacts are sparse over time. Our example is a cattle-trade network, where the movements of animals between farms in Germany are recorded on a daily basis. Next, we isolate the trade within each federal state of Germany as visualized in Fig. 7 and restrict trade to the GSCC of the underlying aggregated graph. Disregarding the smallest networks (those with less than 27 nodes), we thus obtain 12 time-varying graphs with sizes varying from 254 to 27 863 nodes and highly heterogeneous topological and temporal features (see Appendix B). As in the previous section, we assume that premises can be either susceptible, infected, or recovered, and trade events facilitate the transmission of a disease. Unlike before, however, we take into account the number of traded animals during each transaction, i.e., the weight w k→l of a (temporal) contact from node k to l. To this end, we modify the infection propagator in Eq. (15) and replace β by 1 − ð1 − βÞ w k→l (see Appendix A for more information). In a potential outbreak, we assume that an infected node is detected with a constant probability μ each day, after which it would be isolated and thus removed from the network. As a consequence, highly infectious diseases such as FMD can be modeled as SIRtype epidemics [60] . In Fig. 8 , we compare the critical infection probability similar to Fig. 6 for six selected federal states with different transition characteristics. The critical value derived from MC simulations varies between β crit ¼ 0.018 [Bavaria (BY)] and β crit ¼ 1.0 [Saxony (SN)]. The latter indicates that outbreaks remain localized for every choice of β due to sparse intrastate trade. As a potential application to public health institutions, we present in Fig. 9 (a) the spatial variation of the epidemic risk in terms of β crit . The quantitative comparison in Fig. 9 (b) demonstrates that spectral methods provide a lower bound with a varying degree of accuracy depending on the network details. Despite their heterogeneity in size and activity, we find for all networks that the CB model outperforms the IB approach. The detailed results for all states as well as a similar analysis for μ −1 ¼ 120 are available in Appendix C. In this paper, we have presented the CB model for epidemic SIR spreading on temporal networks as a conceptually similar framework to the widely used IB approach. Derived from the message-passing framework [23, 51] , it inherits its accuracy on loop-free topologies and improves analytic estimations with respect to the IB approach for arbitrary time-evolving graphs. Moreover, the focus on edge-based quantities that are updated in discrete time steps allows a seamless integration of temporal interactions. Structurally similar to the node-centric IB model, the proposed CB approach poses a low conceptual barrier and admits application on large graphs. Importantly, the accuracy of the CB model improves existing approximations of the epidemic threshold, which is a crucial risk measure for public health institutions. To this end, we have studied the largest eigenvalue of the infection propagator matrix, which determines the disease propagation in the low prevalence limit and takes into account the full temporal and topological information up to the observation time. The largest eigenvalue can be easily found through repeated matrix multiplications, i.e., the so-called power method. Without relying on extensive MC simulations and a subsequent parameter fit, the critical value can thus be estimated with efficient, vectorizable tools from linear algebra that are available for most high-level programming languages. In the application section, we focused first on a social contact graph that can be used to analyze the propagation of airborne diseases as well as the spread of information. Our comparison between MC simulations and analytic estimations from the CB and IB models followed a bottom-up approach: We looked at (i) epidemic trajectories of individual nodes, (ii) averaged trajectories given the same outbreak location, and (iii) the final outbreak size for a range of infection probabilities and with random initial condition. In all cases, the CB model provides a closer upper bound to MC simulations than the widely used IB model. All results based on the conference data set can be reproduced using the PYTHON code provided in Ref. [57] . Table I ). (b) Critical infection probability β crit in logarithmic scale, sorted from high (left) to low risk (right). Results from MC simulations, and the CB and IB models are presented as groups of blue, red, and green bars, respectively. Disclaimer: A realistic vulnerability analysis requires, for instance, heterogeneous recovery probabilities and complex countermeasures (see Appendix B for details). As a particularly important application, we then compared analytic estimations of the critical infection probability with extensive MC simulations. To this end, we included a case study of livestock trade within 12 federal states in Germany with highly heterogeneous characteristics in terms of size, density, and temporal activity. Consistently, we found that the CB model improves the previously proposed lower bound at a low conceptual and computational cost. Many excellent results have already been derived within the IB framework for empirical networks and in the context of random graphs (see Ref. [69] for a recent review) that can further improve the CB model. We therefore expect that the conceptual simplicity of the CB framework allows us to integrate features such as non-Markovianity [22] , stochastic effects [70] , and estimations of uncertainty [71] that are important to realistic disease models on temporal networks. Also, first steps towards higher-order models that go beyond the tree-graph assumption have been proposed in the context of percolation theory [72] and diffusive transport [73, 74] , and we expect these improvements to be applicable to the CB model as well. The data are available on [75], and using the source code in Ref. [57] , results of this paper can be easily reproduced. In order to improve the predictive power of a network model, it is often required to take into account additional information. In the main article, we focus on the temporal dimension. However, another important piece of information is the weight of a contact. The interpretation of weight can range from passenger numbers in the global air traffic network to the impedance in a network of electric components. The distribution of weights in static as well as temporal empirical networks often shows a broad tail [76] , and as such, the averaged edge weight can become meaningless because of large fluctuations. It is therefore often required to account for heterogeneous edge weights explicitly in epidemiological models. However, depending on the interpretation, weights may enter the model in different ways. Typically, a time-dependent edge weight w k→l ðtÞ is considered similar to the conductivity between two nodes k and l in an electric circuit. Translated to an epidemiological context, we would thus scale the infection probability linearly; i.e., βa k→l ðtÞ becomes βw k→l ðtÞ in a weighted network. Another approach, popular in the context of random walks and disease spreading, is to interpret an integer-valued weight w k→l ðtÞ as a number of parallel and unweighted edges that connect k with l [7, 46] . From an epidemiological viewpoint, this idea would translate to w k→l independent attempts to transmit the disease at time t. Here, the infection probability βa k→l ðtÞ becomes 1− ð1 − βÞ w k→l ðtÞ in the weighted case. In the main article, we apply the latter interpretation to calculate the epidemic risk in the context of livestock trade (see Fig. 9 ). Here, weights correspond to the number of animals traded, each of which can infect the target population independently. For small probabilities β ≪ 1, the adjusted infection probability simplifies to 1 − ð1 − βÞ w k→l ðtÞ ≈ βw k→l ðtÞ. A second source of heterogeneity that is commonly considered includes heterogeneous infection and recovery probabilities, denoted as β k→l and μ k , respectively. With these modifications, the dynamic equations (5)-(7) from the main text translate to θ k→l ðt þ 1Þ ¼ θ k→l ðtÞ − Ψ k→l ðtÞI k→l ðtÞ; ðA1aÞ Here, Ψ k→l ðtÞ denotes the probability that k infects l at time t given that the former is infected and has not yet transmitted the disease. For weighted networks, we can choose Ψ k→l ðtÞ ¼ β k→l w k→l ðtÞ or Ψ k→l ðtÞ ¼ 1 − ð1 − β k→l Þ w k→lðtÞ as discussed above. The linearization of Eqs. (A1a)-(A1c) around the disease-free state leads to Here, the circle denotes the elementwise product. Moreover, the L-dimensional vectors μ and ΨðtÞ integrate the nodeand edge-dependent values μ k and Ψ k→l ðtÞ, respectively. We also generalize the temporal nonbacktracking matrix B β ðtÞ from Eq. (13) to the weighted one: The largest eigenvalue ρ of the infection propagator determines the asymptotic stability for small perturbations around the disease-free state. Accounting for heterogeneity in β, μ, and contact weights, the criticality condition Eq. (16) from the main text reads Assuming β ≡ β k→l for all edges k → l, we can determine from Eq. (A4) the critical (homogeneous) infection probability β crit given a weighted, temporal network with heterogeneous recovery probabilities μ k . Similarly, one can assume μ ≡ μ k for all nodes k and thus derive the critical (homogeneous) value μ crit with heterogeneity in the infection probability β k→l . The system of traceability of cattle in the EU requires that each animal is identified with ear tags and that each movement, birth, or death event has to be reported within 7 days of the event to the national livestock database. We consider an excerpt of the national German livestock database HI-Tier [59] for the year 2010. The database is administered by the Bavarian State Ministry for Agriculture and Forestry on behalf of the German Federal States. It records 3.2 million animal movements with a total of 13.4 million traded animals between 183 454 premises, such as farms, pastures, slaughterhouses, and traders within the observation window. The location of each animal holding was provided at the resolution of the municipality. We consider each trade event between two premises a temporal contact, and we identify an edge if at least one contact has been recorded. The distribution of edges is highly heterogeneous in terms of geography, degree, and weight. In Fig. 10(a) , we observe clusters of trade activity mostly within and between North Rhine-Westphalia (NW), Lower Saxony (NI), Baden-Württemberg (BW), and Bavaria (BY). The number of trading partners, i.e., the node degree, is broadly distributed as demonstrated in Fig. 10(b) . Here, we differentiate between in, out, and total degree. Similarly, we find a broad distribution of edge weights in Fig. 10(c) , i.e., the number of traded animals along a given edge. The geographic distribution of nodes in Fig. 11 (a) shows dense regions in the northwest and southeast including the above-mentioned federal states NW, NI, BW, and BY. Here, we also find the largest premises in terms of total traded animals: In Fig. 11(c) , color and size indicate the node strength, i.e., the aggregated trade volume. The heterogeneous distribution of strength also becomes apparent in Fig. 11(d) where in, out, and total strength are analyzed separately. Finally, we observe in Fig. 11 (b) the net flux, i.e., the difference between in-and out-directed trade volume. We display only nodes with at least 500 traded animals in Figs. 11(b) and 11(c) . From a temporal perspective, we find that trade fluctuates between 10 2 and 10 4 active nodes, i.e., farms with at least one trade event on a given day, whereas minima appear regularly on the weekends [see Fig. 12(a) ]. The weekly pattern is also apparent in the interactivation time distribution, i.e., the time interval between two successive trade events for a given node (see Fig. 12 ). Here, we find a broad distribution of activity with peaks around 7, 14, and 21 days. The geographic risk analysis in Sec. V B 1 requires us to separate the network into subgraphs that correspond to the intrastate trade (see Fig. 7 ). The largest eigenvalue of the corresponding infection propagator allows us to evaluate the outbreak risk within a federal state due to the local movement of infected animals. In Table I , we list the names of all 16 federal states of Germany together with the corresponding ISO abbreviation and basic statistics: the number of nodes, (static) edges, and (temporal) contacts in the GSCC. The city states Berlin, Hamburg, and Bremen as well as Saarland, which is a particularly small state in terms of nodes, are marked with an asterisk and are not considered for risk analysis in Fig. 9 . Separating the trade network into subgraphs as visualized in Fig. 7(b) inevitably reduces the outbreak risk as the neglected cross-border edges would otherwise facilitate the disease transmission. In Fig. 13(a) , we find that a considerable fraction of trade is directed across federal states and has thus been removed. This case applies, in particular, to the federal states NI, NW, and BW. Similarly, we find that the ratio between intrastate and in-directed trade lies between 0.6 (NW) and 0.9 (BY). Thus, we conclude that a considerable fraction of trade across borders is being neglected in the geographic risk analysis in Fig. 9 . It is also important to stress that we use the same parameter μ across all federal states and thus assume a uniform detection probability. In reality, federal states with a large number of premises tend to enforce stricter hygiene and intervention standards so that the actual epidemic risk for states such as BY and NW is much lower. A realistic evaluation for public health must therefore include heterogeneous recovery (detection) probabilities on the level of states or individual nodes as discussed in Ref. [77] and Appendix A, as well as a complex disease response that Table I includes trade restrictions, increased awareness, and higher biosecurity. From the detailed, node-level infection trajectory, we can estimate the infection arrival time from a given outbreak location to all remaining nodes. For that purpose, we extend the contact sequence periodically in time until the infection probabilities are negligible. Then, we derive the infection arrival time to a single node from the corresponding cumulative infection probability (see Fig. 2 ) as follows: (i) The discrete derivative of the cumulative infection probability gives the probability distribution to contract the infection at a given time step. (ii) The expectation value of the probability distribution gives the mean infection arrival time at a single node, corresponding to a scatter point in Fig. 14(a) . Here, we compare the expected values from MC simulations with the estimated infection arrival times given by the CB and IB models, respectively. In a perfect prediction, the scattered values would lie on the diagonal, but, as the contact network is far from a treelike structure, the models estimate infection arrival times smaller than the observed values. The comparison between the IB and CB frameworks in Fig. 14(b) shows a considerably smaller relative deviation of CB estimations from the corresponding MC simulations for the given set of disease parameters and outbreak locations. Another application focuses on the vulnerability of nodes with respect to a given outbreak location. Again, we assume an infinite time horizon and compare the cumulative probability that a node has been infected in the limit t → ∞. As before, we find a good correlation between simulations and the estimated vulnerability in Fig. 15 , whereas the CB model consistently outperforms the IB approach and overestimates the expected values surprisingly little given that the underlying aggregated network is fairly dense (the average degree is hki ≈ 19) and far from being treelike. For some applications, we may be interested in the trajectory of a global epidemic, averaged over outbreak locations. A sufficiently large number of initially infected nodes would then avoid complications with the early outbreak phase [39, 78] . In this case, we adjust the MC simulations such that every node is infected independently with a given probability 1 − z l ¼ 1 − z, ∀ l ∈ N at t ¼ 0. As for the analytic approach, we only need to set a corresponding homogeneous initial condition, and thus the computational complexity remains the same as in the previous case of one initially infected node. In Fig. 16 (left column), we observe a narrow, timedependent distribution of cumulatively infected nodes around the mean value for three different infection probabilities. Without applying any additional threshold, we find a close agreement between the averaged trajectory, and the CB and IB models in all cases. In contrast to Fig. 3 of the main text, we observe in Fig. 16 (right column) only one peak in the distribution due to the large number of initially infected nodes. One potential application is to calculate the vulnerability of a node as discussed in Ref. [6] . Here, the vulnerability is defined as the probability that a given node is eventually infected by a disease that started somewhere in the network. The value can be used to rank nodes in order to prioritize surveillance or vaccination measures to the nodes that are most likely to contract the disease when resources are limited. In Fig. 17 , every curve represents the vulnerability FIG. 14. (a) Comparison between simulated and estimated mean infection arrival times. We extend the data set periodically in time until the outbreak dies out. The discrete derivative of the cumulative infection probabilities (see Fig. 2 ) yields the infection arrival probabilities of which we take the average value for every node. Results according to the CB and IB models are visualized as red circles and green crosses, respectively. The epidemic starts from the same outbreak origin and disease parameters as in Fig. 1 . 15. (a) Comparison between simulated and estimated vulnerability. We compute the cumulative infection probability in the limit t → ∞, also denoted as vulnerability. The comparison with CB and IB estimations is visualized by red circles and green crosses, respectively. Each value corresponds to the vulnerability of a node given the same outbreak location and disease parameters as in Fig. 1. (b) Relative deviation of the estimated values with respect to MC simulations. The numerical values are averaged over 10 5 realizations. of one node as a function of the infection probability β. These results are derived from the CB model, given an initial infection probability of 1 − z ¼ 0.2. The individual colors correspond to the degree of the node in the underlying time-aggregated graph and serve as a guide to the eye. Interestingly, we find that the ranking, as estimated by the CB model, may change with increasing infection probabilities β as can be seen from the highlighted curve in Fig. 17 . This effect has been observed earlier in the context of static networks [6] and indicates that network properties alone are often not sufficient to rank nodes as they do not take into account details of the dynamic system. The analysis of the conference data set in the main text was limited to a single value of the recovery probability with μ ¼ 2.85 × 10 −4 . This choice corresponds to an expected infectious period of about 19.5 h. In addition to the analysis of the main text, we present in Fig. 18 similar results for different values of μ. The left, middle, and right columns correspond to Figs. 4, 6(a), and 6(b), respectively. In all cases, the CB model gives a closer bound to MC simulations as compared to the IB approach. Next, we complete the analysis in Fig. 8 of the main text on the epidemic threshold for the example of the animal trade network. Again, we assume μ ¼ 1=28 and provide in Fig. 19 a detailed analysis of all federal states, excluding the city states Berlin, Hamburg, and Bremen. The left and middle columns in Fig. 19 provide a similar analysis to Figs. 6(a) and 6(b) of the main text, respectively. In other words, we present the distribution of outbreak sizes for different values of the infection probability β (left column) from which we derive the coefficient of variation c v (blue line, middle column). The right column presents values of c v that are close to the peak and a quadratic fit (green line, right column) that determines the numerical estimation of the critical infection probability (blue vertical line). This value can be compared to spectral estimations from the mean-field models. In agreement with previous results, we find that the criticality condition in Eq. (16) of the CB model (see main text) improves previous results of the IB approach. FIG. 16 . Cumulative infection probability with a large fraction of initially infected nodes for three different values of β and the same outbreak location as in Fig. 3(a) . Left column [panels (a1)-(c1)]: Time-evolving distribution (linear scale) of cumulatively infected individuals for infection probabilities β ¼ 0.01, 0.02, and 0.1, respectively. We average over outbreak locations with 20% of the network initially infected at random. The mean outbreak size (blue line), averaged over 10 5 realizations and with a standard deviation below 10 −4 , can thus be compared to the CB and IB models (red dashed and green dotted lines, respectively) with no threshold applied. Right column [panels (a2)-(c2)]: Final distribution (logarithmic scale) together with the averaged values. FIG. 17. Vulnerability as a function of the infection probability β estimated from the CB model. Each curve represents the vulnerability of a node, i.e., the probability to contract the infection from a set of randomly chosen outbreak locations. Here, we estimate the vulnerability according to the CB model. Starting from an initial infection probability of 1 − z ¼ 0.2, we propagate the infection over time until convergence. We stop when the largest increase in vulnerability after 24 h falls below 10 −3 . The colors indicate the degree of each node in the underlying time-aggregated graph. Moreover, the behavior of one selected node is highlighted. Finally, we provide an additional analysis of the epidemic threshold for the cattle-trade data with μ ¼ 1=120. Results in Fig. 20 are akin to our previous analysis in Fig. 19 , except for SL. Here, the spectral condition in Eq. (16) of the CB model predicts that every outbreak remains localized, i.e., β CB crit ¼ 1, whereas MC simulations suggest a transition to global epidemics, hence β MC crit < 1. We attribute the inconsistency to the small size of the network (26 nodes). The spectral approach implicitly assumes an infinitely large network, which is clearly violated in this case. We summarize the results for μ ¼ 1=120 in Fig. 21 , akin to Fig. 9 of the main text. The risk map in Fig. 21(a) visualizes the spatial variability of the outbreak risk, and each group of blue, red, and green bars in Fig. 21 (b) provides a quantitative comparison between MC results, the CB model, and the IB approach, respectively. The Hidden Geometry of Complex, Network-Driven Contagion Phenomena Effective Distances for Epidemics Spreading on Complex Networks Locating the Source of Large-Scale Outbreaks of Foodborne Disease Bayesian Inference of Epidemics on Networks via Belief Propagation Effective Approach to Epidemic Containment Using Link Equations in Complex Networks Assessing Node Risk and Vulnerability in Epidemics on Networks Analytical Computation of the Epidemic Threshold on Temporal Networks Analytical Computation of the Epidemic Threshold on Temporal Networks Containing Epidemic Outbreaks by Message-Passing Techniques On the Definition and the Computation of the Basic Reproduction Ratio R 0 in Models for Infectious Diseases in Heterogeneous Populations Epidemic Spreading in Scale-Free Networks Spread of Epidemic Disease on Networks Epidemic Thresholds in Real Networks Spread of Infectious Disease through Clustered Populations Epidemic Thresholds in Dynamic Contact Networks Percolation on Sparse Networks A Contribution to the Mathematical Theory of Epidemics The Mathematical Theory of Infectious Diseases and Its Applications Exact Epidemic Models on Graphs Using Graph-Automorphism Driven Lumping Virus Spread in Networks Generalization of Pairwise Models to Non-Markovian Epidemics on Networks Mean-Field Models for Non-Markovian Epidemics on Networks Message Passing Approach for General Epidemic Models Oscillations in SIRS Model with Distributed Delays Non-Markovian Infection Spread Dramatically Alters the Susceptible-Infected-Susceptible Epidemic Threshold in Networks Infection Dynamics on Scale-Free Networks Networks and Epidemic Models The Structure and Dynamics of Networks Dynamical Systems on Networks Measuring Large-Scale Social Networks with High Resolution Fundamental Structures of Dynamic Social Networks Temporal Networks Spatial variation of epidemic risk for μ ¼ 1=120. The critical infection probability β crit determines the color of the federal states (see Fig. 20 for details). The visualization is akin to Fig Time-Varying Graphs and Dynamic Networks Information Dynamics Shape the Sexual Networks of Internet-Mediated Prostitution Inferring Friendship Network Structure by Using Mobile Phone Data Diffusion on Networked Systems Is a Question of Time or Structure Epidemic Dynamics on an Adaptive Network Edge-Based Compartmental Modelling for Infectious Disease Spread Infections on Temporal Networks-A Matrix-Based Approach Livestock Disease Spread through Animal Movements: A Static and Temporal Network Analysis of Pig Trade in Germany Network-Based Assessment of the Vulnerability of Italian Regions to Bovine Brucellosis Epidemic Spreading in Real Networks: An Eigenvalue Viewpoint Individual-Based Approach to Epidemic Processes on Arbitrary Dynamic Contact Networks The Effect of Network Topology on the Spread of Epidemics Discrete-Time Markov Chain Approach to Contact-Based Disease Spreading in Complex Networks An Individual-Based Approach to SIR Epidemics in Contact Networks Predicting the Epidemic Threshold of the Susceptible-Infected-Recovered Model Epidemic Processes in Complex Networks Competition-Induced Criticality in a Model of Meme Popularity Inferring the Origin of an Epidemic with a Dynamic Message-Passing Algorithm Spread of Epidemics in Time-Dependent Network Virus Propagation on Time-Varying Networks: Theory and Immunization Algorithms Infection Propagator Approach to Compute Epidemic Thresholds on Temporal Networks: Impact of Immunity and of Limited Temporal Resolution GRAZE: A Generic and Standardized PYTHON Package for the Analysis of Livestock Trade Data Temporal Interactions Facilitate Endemicity in the Susceptible-Infected-Susceptible Epidemic Model Source Code of "Contact-Based Model for Epidemic Spreading on Temporal Networks What's in a Crowd? Analysis of Faceto-Face Behavioral Networks The HIT database (Herkunftssicherungs-und Informationssystem für Tiere) is administered by the Bavarian State Ministry for Agriculture and Forestry on behalf of the German federal states Modeling Infectious Diseases in Humans and Animals Solving the Dynamic Correlation Problem of the Susceptible-Infected-Susceptible Model on Networks Message-Passing Approach for Recurrent-State Epidemic Models on Networks Exact Deterministic Representation of Markovian SIR Epidemics on Networks with and without Loops Dynamic Cavity Method and Problems on Graphs, Thesis The Unreasonable Effectiveness of Tree-Based Theory for Networks with Clustering Introduction to Percolation Theory Predicting Epidemic Risk from Past Temporal Contact Data Mathematics of Epidemics on Networks: From Exact to Approximate Models Temporal Percolation of the Susceptible Network in an Epidemic Spreading The Relationships between Message Passing, Pairwise, Kermack-McKendrick and Stochastic SIR Epidemic Models Beyond the Locally Treelike Approximation for Percolation on Real Networks Causality-Driven Slow-Down and Speed-Up of Diffusion in Non-Markovian Temporal Networks From Networks to Optimal Higher-Order Models of Complex Systems The Architecture of Complex Weighted Networks Disease Persistence on Temporal Contact Networks Accounting for Heterogeneous Infectious Periods Epidemics on Networks with Large Initial Conditions or Changing Structure