Unifying continuous, discrete, and hybrid susceptible-infected-recovered processes on networks

Waiting times between two consecutive infection and recovery events in spreading processes are often assumed to be exponentially distributed, which results in Markovian (i.e., memoryless) continuous spreading dynamics. However, this is not taking into account memory (correlation) effects and discrete interactions that have been identiﬁed as relevant in social, transportation, and disease dynamics. We introduce a framework to model continuous, discrete, and hybrid forms of (non-)Markovian susceptible-infected-recovered (SIR) stochastic processes on networks. The hybrid SIR processes that we study in this paper describe infections as discrete-time Markovian and recovery events as continuous-time non-Markovian processes, which mimic the distribution of cell cycles. Our results suggest that the effective-infection-rate description of epidemic processes fails to uniquely capture the behavior of such hybrid and also general non-Markovian disease dynamics. Providing a unifying description of general Markovian and non-Markovian disease outbreaks, we instead show that the mean transmissibility produces the same phase diagrams independent of the underlying interevent-time distributions.


I. INTRODUCTION
Models of epidemic processes such as the susceptibleinfected-recovered (SIR) model and related models have provided various insights into dynamical and stationary features of disease, opinion, and failure spread in social and technical systems [1][2][3].In the SIR model, infected individuals may transmit a disease to susceptible ones.After a certain period, infected individuals recover and are not part of the diseasetransmission process anymore.The exact time evolution of the continuous-time stochastic SIR spreading process is described by the Chapman-Kolmogorov equation or its differential form (i.e., the master equation).However, exact analytical solutions of the master equation are limited to special cases and therefore different approximations are being used (e.g., deterministic ordinary differential equation models [4][5][6], cavitylike models [7,8], and pairwise approaches [9]).
Gillespie and kinetic Monte Carlo (KMC) approaches [10] provide techniques to generate statistically exact trajectories of a master equation.The assumption underlying KMC simulations of SIR processes is that waiting times between consecutive recovery and infection events are exponentially distributed.However, many natural processes including social dynamics [11,12] exhibit correlation and memory (i.e., non-Markovian) effects [13] and therefore are not described by exponential (i.e., memoryless) waiting-time distributions [14].Recent attempts to simulate general non-Markovian processes led to the development of the non-Markovian Gillespie algorithm (NMGA) [15] and the Laplace Gillespie algorithm (LGA) [16].Both methods are based on a mapping of multiple stochastic processes with general (continuous-time) waitingtime distributions to a modified KMC algorithm.As outlined in Ref. [16], the NMGA is exact only for infinitely many processes and requires one to recalculate all individual rates at every time step.This algorithm has the advantage that it can simulate arbitrary continuous-interevent-time distributions [15].The LGA interprets survival functions of waitingtime distributions as Laplace transforms of underlying eventrate distributions.Although the LGA is exact for arbitrary numbers of processes, it is only applicable to certain waitingtime distributions [16].
In the context of non-Markovian SIR models, there also exist event-driven and directed-percolation-based approaches that are directly applicable to these types of processes [17,18].In this paper we use an approach similar to the directedpercolation method of Refs.[17,18] and map non-Markovian SIR processes with general waiting-time distributions to a shortest-path problem [19,20] in a weighted spreading network.We refer to this method as the shortest-path kinetic Monte Carlo (SPKMC).In contrast to the NMGA, LGA, and aforementioned approximation techniques [4][5][6][7][8][9], our SPKMC framework allows us to produce exact stochastic trajectories of SIR processes for general continuous-and discrete-waitingtime distributions on any network.
We use our framework to study continuous, discrete, and hybrid SIR processes.Such hybrid formulations of spreading dynamics can be useful to account for events that are generated according to a sequence of fixed-duration schedules (e.g., meetings [21], patient care [22,23], and transportation schedules [24]) and finite transfer times of information in communication networks [25].
Various processes that affect human interaction and information exchange in communication networks can be modeled as discrete events.Possible examples include sequences of fixed-duration schedules of meetings (multiples of hours and minutes) [21][22][23] and transportation schedules [24] that can affect disease transmission on a metapopulation level.For communication networks, finite transfer times and synced digital events (e.g., release of computer viruses) are also examples of discrete processes [25].
In addition, hybrid SIR models can also account for latency periods that are typically modeled by introducing an additional "exposed" compartment [5] and a corresponding latent period, which mimics the observed mean incubation time (e.g., 8-14 days for measles [26]).As an alternative to discrete [27][28][29] and continuous [5] epidemic models with exposed compartments, our approach can account for discrete-time delays by directly modifying transmission rates.
Our results suggest that hybrid and also general non-Markovian disease outbreaks cannot be uniquely captured by effective infection rates.However, by mapping hybrid SIR processes to bond percolation [30,31], we show that the corresponding mean transmissibilities produce the same phase diagrams independent of the underlying infection-and recovery-time distributions, thus providing a unifying description of general Markovian and non-Markovian SIR processes.The framework we propose assumes no specific form of waiting-time distributions [29,32] and allows us to simulate and analytically describe general discrete, continuous, and hybrid variants of Markovian and non-Markovian SIR dynamics.

II. SHORTEST-PATH KINETIC MONTE CARLO
Before focusing on the simulation of hybrid SIR dynamics, we introduce the necessary mathematical toolbox that allows us to map general waiting times to a shortest-path problem in an underlying spreading network.Let φ(τ ) and ψ (τ ) be the probability density function (PDF) or probability mass function of recovery and infection times.In continuous time, the probability that a recovery (infection) event occurs in the interval (τ, τ + dτ ) is φ(τ )dτ [ψ (τ )dτ ].The discretetime analogs φ(τ ) and ψ (τ ) are the recovery and infection probabilities after τ steps, where • denotes the floor function.We denote the cumulative distribution functions (CDFs) of φ(τ ) and ψ (τ ) by (τ ) and (τ ).The function (τ ) [ (τ )] is the probability that a recovery (infection) event occurred in [0, τ ].In the case of Poissonian SIR dynamics, waiting-time distributions are described by the PDFs φ(τ ) = γ e −γ τ and ψ (τ ) = βe −βτ and CDFs (τ ) = 1 − e −γ τ and (τ ) = 1 − e −βτ , where γ and β are the corresponding recovery and infection rates.Note that we use prefixes such as Erlang-geometric to indicate the recovery-and infection-time distributions [φ(τ ) and ψ (τ )] of the corresponding hybrid SIR process.
For given distributions φ(τ ) and ψ (τ ), we consider M realizations of SIR dynamics to correspond to an ensemble of M directed spreading networks {G k (V, E )} k∈{1,...,M} , where V and E denote the sets of nodes and edges.Each network G k (V, E ) is initialized as follows.For each node s in G k (V, E ) we generate a random number x ∼ U (0, 1) and use an inverse transform sampling of (τ ) to determine the recovery time of node s using −1 (x).For each node t that is adjacent to s we generate another random number y ∼ U (0, 1) and determine the infection time −1 (y).We now use −1 (x) and −1 (y) to determine edge weights [17] We set ρ st to −1 (y) (i.e., the disease transmission time from node s to t) if infection occurs before recovery and ρ st = ∞ otherwise.Note that the interaction terms can also be general we show an illustration of the weight initialization procedure for a network that consists of four nodes.We again note that networks G k (V, E ) are directed (i.e., weights ρ st may be different from ρ ts ).In the case of Poissonian dynamics, we obtain ( In addition to edge weights, we also keep track of node weights τ i = −1 (x) to describe the evolution of SIR dynamics in a network.
After having identified all weights, we obtain one realization G k (V, E ) of the spreading network.In the next step we take G k (V, E ) and infect one uniformly at random selected node [see Fig. 1(b)].All nodes that are connected to the initially infected node through paths of finite length are also infected and recover in the limit of τ → ∞ [see Fig. 1(c)].The shortest-path length between an infected source and its target node is the corresponding disease transmission time.This formulation of disease transmission can be viewed as a leastaction principle for KMC.If all paths that connect two nodes are infinite, we know that one or multiple recovered (or quarantined or removed) nodes hinder disease transmission (see node 4 in Fig. 1).To describe SIR dynamics with n initially infected nodes, we use I j k to denote the set of infected nodes that result from an initial infection of node j in G k (E , V ).The corresponding set of all infected nodes that result from multiple spreading seeds in G k (V, E ) is I k = n j=1 I j k .The SPKMC framework also allows us to monitor the infection and recovery times of individual nodes.In Appendixes A and B we outline how the evolution of susceptible, infected, and recovered nodes can be reconstructed from shortest paths and describe the possibility to account for quarantine protocols in SPKMC simulations.The stationary fraction of recovered nodes in network is the number of nodes, and the corresponding ensemble average over For each network realization G k (E , V ) we have to identify all shortest paths, an operation with a runtime of order O(|E | + |V | log |V |) when using optimized data structures such as Fibonacci heaps [33].For n initially infected nodes, we can determine the number of susceptible S(t ), infected I (t ), and recovered R(t ) nodes at time t (see Appendix A) after running Dijkstra's algorithm n times.Typically, the number of initially infected nodes n is small and thus the runtime of our algorithm is still of order Reference [16] discusses the runtime complexity of the NMGA and LGA per generated event.Since our SPKMC framework can simulate SIR dynamics with one run of Dijkstra's algorithm, the computational complexity of our framework does not depend on the number of time steps that one wants to simulate.
In Appendix C we consider Poissonian dynamics [see Eq. ( 2)] and show that our SPKMC simulations of stationary and dynamical SIR features agree well with corresponding KMC simulations.An advantage of our proposed shortest-path SIR simulation method is the possibility to simulate general Markovian and non-Markovian dynamics with continuous-and discrete-waiting-time distributions.

III. HYBRID CONTINUOUS-DISCRETE SIR DYNAMICS
To describe latency periods in infection processes (i.e., no infection occurs during a certain time window), we apply our simulation framework to hybrid Poisson-geometric SIR dynamics with discrete infection events that are distributed according to a geometric density function where q is the probability that an infection event occurs within a time step of 1 and δ(•) is the Dirac delta function.We show in Appendix D that the master equation of hybrid Poissongeometric SIR dynamics is no longer time homogeneous.
In Appendix E we also consider an alternative definition of the geometric distribution ψ G 2 that takes on finite values for all non-negative integers.The geometric distribution is the discrete memoryless counterpart of exponential distributions.
We can now use our simulation framework to study disease outbreak characteristics of such hybrid SIR processes.As for many epidemic processes [5], we characterize disease dynamics in terms of the effective infection rate λ = τ φ / τ ψ , where τ ψ and τ φ are the mean times to infection and recovery, respectively.For a fully Poissonian SIR process, the effective infection rate is λ PP = β/γ and invariant upon rescaling of infection and recovery rates by a constant factor [1].That is, the corresponding fraction of recovered [see Eq. ( 3)] only depends on the effective infection rate λ PP [see Fig. 2(a)].By analogy, we now use λ PG 1 = 1/γ τ ψ to denote the effective infection rate for Poissonian-geometric dynamics with infection-time PDF ψ G 1 and (5) However, unlike in fully Poissonian SIR dynamics, we cannot uniquely capture the corresponding phase space by the effective infection rate λ PG 1 [see Fig. 2(b)].That is, we observe different fractions of recovered r for the same value of λ PG 1 .
To better understand the phase space of hybrid SIR processes, we proceed with a mapping to bond percolation.

IV. MAPPING HYBRID SIR DYNAMICS TO BOND PERCOLATION
We now analytically characterize the hybrid SIR disease prevalence in terms of the mean transmissibility T that describes the probability of an infection to be transmitted from an infected to an adjacent susceptible node [31]: In networks with no degree correlations, the critical transmissibility above which an SIR epidemic spreads through a finite fraction of the system is given by the bond percolation threshold [1,31] where k and k 2 denote the first and second moments of the degree distribution P k .For Poisson-geometric SIR dynamics, the mean transmissibility is In Appendix E we provide further details about the derivation of Eq. ( 8) and compare the fully Poissonian and Poissongeometric cases.According to Eqs. ( 7) and ( 8), we find the phase separation line q c = p c (e γ − 1)/(1 − p c ), which separates the phases with and without disease outbreaks (see Fig. 3).For γ = 1 and 0.1, we obtain the critical effective infection rates λ c PG 1 ≈ 0.57 and 0.35, respectively.These values agree well with the numerical data of Fig. 2(b).
Note that TPG 1 (γ , q) cannot be parametrized in terms of an effective infection rate λ PG 1 (see Fig. 3) whereas for fully Poissonian dynamics the mean transmissibility T PP = λ PP /(1 + λ PP ) only depends on the effective infection rate λ PP = β/γ (see Appendix C and Ref. [31]).However, if γ is small, we find that Thus, for sufficiently small values of q and γ (i.e., long mean infection and recovery times), the mean transmissibility TPG 1 (γ , q) only depends on the effective infection rate λ PG 1 .
In Appendix E we show that this is also the case for the alternative formulation of Poisson-geometric SIR dynamics with ψ G 2 .A graphical interpretation of this result is that the phase separation lines in Fig. 3 merge as τ φ = γ −1 and τ ψ tend to infinity.
To determine the relative size of the epidemic S ( T ) as a function of the mean transmissibility T , we use a generatingfunction approach [31] and first consider an uncorrelated network for which the conditional probability P(k|k ) = kP k / k does not depend on k .This approach is based on two generating functions G 0 (x; T ) and G 1 (x; T ).The former is the generating function of the distribution of occupied edges belonging to a certain node, as a function of T [31]: The distribution of occupied edges leaving a node at which we arrived by following a randomly selected edge is generated by [31] Note that we use G 0 (x; T ) to indicate a derivative of G 0 (x; T ) with respect to x.To determine S ( T ) [see Eq. ( 13)], we solve the self-consistency equation and compute where u is the probability that the node at the end of a randomly selected edge does not lead to a giant macroscopic component.In Appendix F we generalize Eqs.(11) and (13) to account for correlation effects between nearest neighbors.Note that the generating-function formalism is useful for cases when the exact network is unknown but its degree distribution is known.For details on limitations of the described bondpercolation mapping, see Appendix H and Refs.[8,34,35].
As in Fig. 2, we now consider a random-regular graph with degree k = 5.The degree distribution is P k = δ k5 , where the Kronecker delta is δ kk = 1 if k = k and zero otherwise.In Fig. 2 we show the analytical solution of Eqs. ( 13) and (11) for fully Poissonian ( TPP ) and Poisson-geometric ( TPG 1 ) SIR dynamics.For a small number of initially infected nodes, the relative outbreak size S( T ) corresponds to the fraction of recovered nodes r [see Eq. ( 3)].We observe that the bond-percolation descriptions of hybrid and fully geometric SIR [see Fig. 4(a)] outbreaks agree well with simulations.

V. UNIFYING NON-MARKOVIAN SIR PROCESSES
Up to this point we have focused on hybrid SIR processes with variations in the infection-time distributions.To understand the general applicability of our framework, we now consider non-Markovian SIR dynamics with recovery times that are distributed according to the Erlang distribution where n and γ are the so-called shape and rate parameters.The Erlang distribution allows us to account for recovery processes that are not just exponentially distributed but more concentrated within a certain time window.It is the distribution that describes the sum of n independent exponential variables with rate γ .The Erlang distribution has been used as an approximation of cell-cycle time distributions [36] and as such it is a good candidate for disease recovery processes as cells cycles have n stages through which they are progressing (e.g., n = 4 for COVID-19 [37]).In Fig. 4(b) we observe that Erlang-geometric SIR processes also cannot be described by an effective infection rate.The considered examples of non-Markovian SIR processes show that the effective-spreadingrate description cannot uniquely characterize hybrid and general non-Markovian disease outbreaks.Instead, the mean transmissibility T provides a unifying control parameter as we show in Fig. 4(c).In Appendix G we outline that this also holds for other networks including Erdős-Rényi, Barabási-Albert, and various empirical networks.

VI. DISCUSSION AND CONCLUSION
We introduced numerical and analytical frameworks for the study of general (non-)Markovian SIR dynamics on networks.Furthermore, we proposed a hybrid SIR process that models infection and recovery as discrete-time Markovian and continuous-time (non-)Markovian processes, respectively.
The discussed examples of hybrid SIR processes can account for cell-cycle distributions and latency intervals during which no infection events occur.We showed that the effectiveinfection-rate description of Markovian SIR processes [1] cannot uniquely capture non-Markovian epidemic outbreaks.However, our results suggest that the mean transmissibility provides a unifying description of (non-)Markovian SIR processes across a wide range of network structures (see Appendix G) and infection-and recovery-time distributions.These observations are of particular interest for disease control and hint at a redefinition of the epidemic threshold to appropriately account for disease dynamics and network structure [38].Our results also complement an earlier study on non-Markovian SIR dynamics [39], which showed that strong temporal heterogeneity in the contact patterns between individuals may significantly suppress epidemic outbreaks.
Further motivation for the study of discrete interaction processes comes from temporal-network theory, where the majority of temporal interactions is considered to be discrete [40,41].Future studies may extend our work to hybrid models on temporal networks.
Our findings are in accordance with earlier results on non-Markovian susceptible-infected-susceptible (SIS) dynamics [14], where a modified effective infection rate was used to uniquely capture corresponding steady states.A mean-field analysis of SIS dynamics [42] also revealed that there is an equivalence between certain non-Markovian and Markovian SIS processes.Similar concepts may be helpful to better understand similarities between non-Markovian and Markovian SIR dynamics.
To summarize, our work can contribute to more accurate and informative models of spreading processes on networks and metapopulation spreading models [43,44].

APPENDIX A: DYNAMICS RECONSTRUCTION
Based on the SPKMC framework that we outlined in Sec.II, we can also determine the evolution of susceptible, infected, and recovered nodes.We use 1 − S(t ) to denote the mean number of nonsusceptible nodes prior to some time t and compute this quantity from the set of edge-weighted spreading graphs {G 1 , . . ., G M } according to where I is the set of initially infected nodes and d G i (k, j) is the shortest-path length from node k to node j in the weighted spreading network G i .To determine the fraction of infected and recovered nodes at time t, we compute the cardinality of the set of all nodes that are connected with an infected source node through a path of maximum length t.To extract the dynamical behavior of infected nodes from an SPKMC simulation, we need to include node-recovery weights {τ i } in our simulation framework.Similarly to Eq. (A1), we determine the mean number of infected nodes I (t ) according to

APPENDIX B: QUARANTINE DYNAMICS
We now describe the possibility to apply the SPKMC framework of Sec.II to quarantine modeling.Recall that Dijkstra's algorithm [45] is constructing shortest paths via dynamic programming updates.At iteration n = 0, all distances d G i (k, j) [n] = ∞ are set to infinity, except for the source node d G i (k,k) [n] = 0.At iteration n, we update the shortest path with the equation where ρ j,l is the edge weight between nodes j and l.
One possible quarantine strategy would be that a predefined set of nodes gets removed from the network at quarantine time t 1 and brought back at time t 2 .This procedure can be repeated as often as necessary and directly incorporated in our simulation framework.Nodes can only infect others or be infected by surrounding nodes if they are not under quarantine (i.e., not removed).The outlined quarantine protocol can be implemented as follows.For a given source node k, some target node l, and a node j that is under quarantine during t ∈ [t 1 , t 2 ], the shortest-path calculation can be extended by adapting distance updates in Dijkstra's algorithm where ρ j,l is the interevent edge transmission delay between nodes j and l, B = R 0 \ [t 1 , t 2 ], and χ B (x) is the characteristic function of B: According to Eq. (B5), we obtain a finite shortest-path length for at least one node j.We show an SPKMC simulation for quarantine on a random network in Fig. 5.We observe a drop in the number cumulative proportion of infections c(t ) = t 0 i(t )dt as soon as quarantine begins.
The prior quarantine protocol can be generalized as follows.Each node j can have its own quarantine from time t ( j) 2 .Accordingly, we can define 2 ] and χ B ( j) (x), the characteristic function of B ( j) : Now, for a given source node k, some target node l, and a node j that lies on the path between k and l, the shortest-path , the initial fraction of 10 −3 infected nodes spreads through the system.We show the fraction of recovered nodes r as a function of λ PP .The blue solid line is a solution of the SIR percolation problem [31] [see Eqs. ( 12) and ( 13)] and the gray squares and black crosses are KMC and SPKMC simulations, respectively, averaged over 10 3 realizations (error bars are smaller than the markers).calculation in Dijkstra's algorithm can be extended by adapting the dynamic programming update where the characteristic function χ B ( j) of node j prohibits transmission from node j if it is under quarantine and the characteristic function χ B (l ) prohibits transmission to node l if it is under quarantine.

APPENDIX C: COMPARISON OF ALGORITHMS
Here we compare the stationary and transient behaviors of fully Poissonian SIR dynamics that we obtain with SPKMC and KMC simulations [18] (see Fig. 6).We consider a regular random network with degree k = 5 and N = 10 4 nodes.For Poissonian infection and recovery times, the mean transmissibility is [31] where the effective infection rate is λ PP = β/μ.This yields the threshold above which giant (outbreak) components are observable.We show a comparison of the transient behavior of SPKMC and KMC simulations in Fig. 7.

APPENDIX D: MASTER EQUATION FOR HYBRID POISSON-GEOMETRIC SIR DYNAMICS
In this Appendix we formulate the master equation for hybrid Poisson-geometric SIR dynamics.Let us denote the probability of finding a system in configuration σ at time t by P(σ, t ).The probabilities P(σ, t ) fulfill the normalization condition σ P(σ, t ) = 1 and the probabilistic evolution of the system is governed by the master equation The first term on the right-hand side of Eq. (D1) describes the inflow into configuration σ from other configurations σ * with transition rate W (σ * → σ ) and the second term accounts for the corresponding outflow with transition rate W (σ → σ * ).In the case of SIR dynamics, every configuration is an n-dimensional vector σ = (σ 1 , . . ., σ n ) that describes the state of every node with σ i ∈ {S, I, R}.For hybrid Poisson-geometric SIR process, the transition rates are not constant in time anymore.This implies that the stochastic process is still Markovian, but not time homogeneous as some transitions may only occur for integer times t.We factorize the transition rates as where A is the adjacency matrix and w(•) denotes the local transition rate of state σ i conditioned on the neighboring states σ * j at time t.The recovery rates are constant in time: However, disease transmissions only occur when the time t is an integer where 1 N + (t ) denotes the indicator function of the positive natural numbers N + , which is equal to 1 when t belongs to N + and zero otherwise.

APPENDIX E: INFLUENCE OF DISCRETE INFECTION TIMES ON TRANSMISSIBILITY
The definition of the geometric distribution in Sec.III [see Eq. ( 4)] yields the mean transmissibility where φ(τ ) = γ e −γ τ is the exponential recovery-time distribution with recovery rate γ .For small γ , the mean transmissibility is (E2) where λ PG 1 = 1/γ τ ψ = q/γ .Based on the definition of the counting process in the Bernoulli trials, one can also define the geometric distribution where the counting starts at 1.There is no correct way; it depends on the actual definition of the stochastic process and alignment of the discrete and continuous interevent times.Different definitions have different interpretations.For example, the above definition assumes that we allow the transmission to happen with probability q at any non-negative integer number.However, if the counting process is defined over natural numbers (positive integers), it would represent the process with density and mean transmissibility Note that the limiting behavior lim γ →∞ TPG 2 (γ , q) = q is in sharp contrast to the Poissonian case, where lim γ →∞ TPP (λ PP ) = 0. We find for small values of q and γ , ) where λ PG 2 = 1/γ τ ψ = q/(1 − q)γ .We show a comparison of the phase diagrams of hybrid SIR processes with infection-time distributions ψ G 1 and ψ G 2 in Fig. 8. FIG. 8. Influence of geometric distributions on transmissibility.We illustrate the dependence of the mean transmissibilities (a) TPG 1 (γ , q) and (b) TPG 2 (γ , q) on the infection probability q and recovery rate γ .The black solid line separates the phases with and without giant (outbreak) components.

APPENDIX F: CORRELATED NETWORKS
The bond-percolation mapping that we outlined in Sec.IV is applicable to networks with an uncorrelated degree distribution [i.e., P(k|k ) = kP(k)/ k ].For networks with degree correlations, we use the notation [46] P(k|k ) = k P(k, k )/k P(k ) and extend the generating-function approach of Sec.IV according to where This extension accounts for correlations between neighboring nodes and is based on the assumption that the considered network is locally treelike [46].We use k cut to denote the largest degree in the network.The function G 0 (u; T ) is the generating function of the distribution of occupied edges belonging to a certain node.The distribution of occupied edges leaving a node at which we arrived by following a randomly selected edge is generated by G 1 (u; T ) and u k is the probability that a node with degree k at the end of a randomly selected edge is occupied.Furthermore, note that degree-degree correlations become irrelevant for the percolation transition if the spectrum of the branching matrix B k,k = (k − 1)P(k |k) satisfies certain conditions [46].SIR dynamics on a corresponding network.The generatingfunction approach provides a possibility to gain insights into (non-)Markovian disease outbreaks if we only know about the degree distribution of a certain network.In this case, the underlying network structure would be implicitly described by a configuration model, which corresponds to a network reconstruction using the maximum-entropy principle [47] with a certain degree distribution as constraint.We summarize these points in Fig. 9.

APPENDIX G: DIFFERENT NETWORKS
In the main text we outlined that effective infection rates cannot uniquely capture hybrid and general non-Markovian disease outbreaks.However, our results for random-regular networks (see Fig. 4) suggest that the mean transmissibility [see Eq. ( 6)] produces phase diagrams that are independent of underlying infection-and recovery-time distributions.In Fig. 10 we show that this observation can also be made for other synthetic and real-world networks.For Erdős-Rényi networks, we can use the bond-percolation description of disease outbreaks to analytically describe the phase diagram.In the case of Barabási-Albert networks, we use the generatingfunction approach for correlated networks (see Appendix F) to obtain the corresponding analytical description.

APPENDIX H: FURTHER CORRECTIONS TO THE PERCOLATION MAPPING
We utilized the SPKMC framework and Gillespie methods [10,51] to generate exact realizations of hybrid and non-Markovian SIR dynamics and compared them to corresponding analytical (bond-percolation) predictions on synthetic and real-world networks.The discussed mapping to bond percolation can be further enhanced by introducing certain corrections.For example, corrections to the mean transmissibility may result from considering the probability that m out of n edges get activated [17]: (H1) Note that the mean transmissibility of the main text is a special case of T n,m for m = n = 1.In addition, corrections accounting for semidirected spreading involve tracking down the exact direction of activated links [34] using an extended generating-function formalism.

FIG. 1 .
FIG. 1. SIR dynamics on a spreading network.(a) Initialization of edge weights according to Eq. (1) in a spreading network that consists of four nodes.Finite edge weights indicate that disease transmission can occur along the corresponding edge.(b) Node 1 is infected and transmits the disease to susceptible nodes 2 and 3 that are connected via paths of finite length (indicated by black arrows).The times at which nodes 2 and 3 get infected are ρ 12 = 2.6 and ρ 12 + ρ 23 = 7.5.All paths that connect node 1 with node 4 have infinite length, so node 4 cannot be infected.(c) In the long-time limit τ → ∞, all infected nodes recover.

FIG. 2 .
FIG. 2. Fraction of recovered nodes in Poissonian and hybrid Poisson-geometric SIR processes.The fraction of recovered nodes r [see Eq. (3)] is plotted as a function of the effective infection rate for (a) fully Poissonian and (b) Poisson-geometric SIR processes.The effective spreading rates are λ PP = β/γ and λ PG 1 = q/γ .Analytical solutions (blue solid lines) are based on Eqs.(13) and (11) for a random-regular graph with k = 5.Numerical simulations have been performed for N = 10 5 nodes and M = 10 2 realizations.

FIG. 3 .
FIG. 3. Comparison of phase spaces.We show the separation lines between phases with and without disease outbreaks.The black solid line corresponds to the Poissonian SIR dynamics with T PP = λ PP /(1 + λ PP ) and λ PP = β/γ .The light (dark) gray solid line describes the phase separation for the Poissonian-geometric SIR dynamics with ψ G 1 (ψ G 2 ) [see Eqs.(4) and (8) and Appendixes D and E for details].

FIG. 4 .
FIG. 4. Outbreak sizes for Markovian and non-Markovian.The fraction of recovered nodes r [see Eq. (3)] is shown for (a) geometric-geometric and (b) Erlang-geometric SIR processes as a function of the corresponding effective infection rates.(c) Markovian and non-Markovian outbreak sizes collapse onto the same curve when plotted against T .Numerical simulations have been performed for N = 10 5 nodes and M = 10 2 realizations with 100 random initial infections.
Fellowship on "Multispecies interacting stochastic systems in biology" and the Army Research Office (Grant No. W911NF-18-1-0345).N.A.-F.acknowledges financial support from SoBigData++ through Grant Agreement No. 871042.L.B. and N.A.-F.contributed equally to this work.

FIG. 7 .
FIG. 7. Transient behavior in SPKMC and KMC SIR simulations.We show the fraction 1 − s(t ) of nonsusceptible nodes as a function of time for fully Poissonian SIR dynamics (blue solid line, KMC; black crosses, SPKMC).Simulations have been performed on a random-regular network with degree k = 5 and N = 10 4 nodes and the numerical data are based on M = 10 2 samples.We used a recovery rate γ = 0.1 and infection rates (a) β = 0.05 and (b) β = 0.08.
FIG. 9. Schematic diagram of the proposed framework.(a) Example of a contact network (Barabási-Albert) with N = 225 nodes, each new node is connected to two existing nodes, and node size scales with betweenness centrality.The generating-function formalism is shown for a given network degree distribution (b) without and (c) with degree correlations.(d) SPKMC simulations on a given network.(e) (Non-)Markovian temporal interactions for a generalized SIR process.