Temporal Gillespie algorithm: Fast simulation of contagion processes on time-varying networks

I. INTRODUCTION

• Up to ~500 times faster than rejection sampling for the investigated processes and parameter ranges, the temporal Gillespie algorithm was numerically validated against that method.The comparison covers the SIS and SIR models described in the introduction.Source: §I · p. 2
• The algorithm is extended to non-Markovian processes, including SIR models with non-constant recovery rates.The paper notes that contagion can exhibit memory effects, such as recovery rates that change over time.Source: §I · p. 2

A. Stochastic processes on time-varying networks

• The temporal algorithm handles changes in possible transitions as the network changes, requiring only finitely many changes to the transition set in any finite interval.Its validity also assumes independent transition processes, which may still depend on earlier transitions and node states.Source: A. Stochastic processes on… · p. 2
• The examples distinguish spontaneous, contact-dependent, and mixed transitions, whose rates or updates depend on node state, contacts, or both.This distinction lets spontaneous processes be updated after transitions and mixed processes after relevant contacts, rather than at every time-step.Source: A. Stochastic processes on… · p. 3

B. Rejection sampling for Monte Carlo simulations

• Choose a sufficiently small time-step so discretization does not significantly affect the process; for temporal networks, measurement or simulation resolution often determines it.The passage requires λ_mΔt ≪ 1 to avoid discretization errors.Source: B. Rejection sampling for Monte Carlo… · p. 5
• At each time step t = 0, Δt, 2Δt, ..., draw a uniform random number rm on [0, 1) for each possible transition m and trigger it when rm < λmΔt.For finite Δt, the passage says to compare rm to 1−exp(λmΔt) to ensure λm defines a proper transition rate; the procedures are equivalent as Δt → 0.Source: B. Rejection sampling for Monte Carlo… · p. 5
• Because λ_mΔt must be small, the vast majority of rejection-sampling trials are rejected, making the algorithm computationally inefficient.The rejected-trial proportion is a weighted average over 1−λ_mΔt.Source: B. Rejection sampling for Monte Carlo… · p. 5

C. Gillespie algorithm on static networks

• Between events, the possible-transition set and rates can be treated as constant, giving a cumulative transition rate Λ equal to the sum of the rates.This condition underlies the static-network waiting-time calculation.Source: C. Gillespie algorithm on static… · p. 5
• Draw the waiting time to the next transition from an exponential distribution with rate Λ.The passage writes this as τ ∼ Exp(Λ).Source: C. Gillespie algorithm on static… · p. 5
• Choose the next transition m with probability proportional to its rate, π_m = λ_m/Λ.This selection follows the draw of the inter-event waiting time.Source: C. Gillespie algorithm on static… · p. 5

D. Temporal Gillespie algorithm

• Drawing τ′ from Exp(1) and solving L(t; t*) = τ′ locates the next event despite a time-varying set of transitions and rates.The cumulative rate is integrated over time to define the normalized waiting time.Source: D. Temporal Gillespie algorithm · pp. 5–6
• The algorithm is presented with pseudocode and a C++ implementation for contagion processes on time-varying networks.The implementation is identified as applying to contagion processes on time-varying networks.Source: D. Temporal Gillespie algorithm · p. 6
• At each time step, only contact-dependent and mixed processes tied to contacts at that step need updating; after an event, all transition processes are updated.This distinction lets the algorithm avoid updating every process at every time step.Source: D. Temporal Gillespie algorithm · p. 6
• By construction, simulated waiting times follow their exact distributions, whereas rejection-based sampling is exact only as λ_mΔt → 0.The authors term this property “stochastically exact.”Source: D. Temporal Gillespie algorithm · p. 7
• For contagion networks, tracking individual waiting times in the next reaction method scales with nodes and contacts, while temporal networks also require frequent internal-time updates.The authors argue that these updates create substantial overhead between events, favoring a single global normalized waiting time.Source: D. Temporal Gillespie algorithm · p. 7

Numerical validation

• For sufficiently small λ_mΔt, SIR and SIS outcomes from temporal Gillespie and rejection sampling are indistinguishable; outside that regime, rejection sampling may be inaccurate.The temporal Gillespie algorithm is stochastically exact for all λ_mΔt, whereas rejection sampling is expected to be accurate only when λ_mΔt ≪ 1.Source: E. Comparison of Gillespie and… · p. 7, Figure 3 · p. 8

Comparison of simulation speed

• The temporal Gillespie algorithm is up to multiple orders of magnitude faster than rejection sampling, with larger gains on larger and sparser networks.The comparison covers SIR and SIS simulations, including empirical time-varying networks.Source: E. Comparison of Gillespie and… · p. 7, Figure 4 · p. 9, Figure 5 · p. 10
• Speed gains are larger on sparser networks, but when all transitions depend on contacts, temporal Gillespie is only approximately twice as fast as rejection sampling.Finding susceptible-infected contacts to determine possible infections limits temporal Gillespie's performance in that case.Source: E. Comparison of Gillespie and… · pp. 7–8

Expected time complexity of the algorithms

• For sparse networks with small transition rates, rejection sampling is limited by rejected trials, while temporal Gillespie is limited by scanning contacts to update possible transitions.The performance gap decreases as mean instantaneous degree increases.Source: E. Comparison of Gillespie and… · p. 9

Improving performance by removing obsolete contacts

• Removing obsolete contacts in empirical SIR simulations pays off when the process's slowest timescale exceeds data length, but can slow simulations when the process is faster.For fast processes the data is repeated few times, so contact-removal overhead can outweigh its savings; for slow processes repeated scans are reduced.Source: E. Comparison of Gillespie and… · pp. 9–10, Figure 5 · p. 10

Slow network dynamics

• When network dynamics are much slower than the spreading process, temporal Gillespie behaves like static Gillespie between network changes and accounts for changes when they occur.In this regime, rejection sampling requires a simulation time-step much smaller than the network-change time-step.Source: E. Comparison of Gillespie and… · p. 10

F. Non-Markovian processes

• Non-Markovian transition rates depend on process history, and independent transition processes determine the waiting-time survival probability.The histories record information relevant to each process up to the current time.Source: F. Non-Markovian processes · p. 11
• The next transition time requires solving an implicit equation when rates vary continuously; treating the cumulative rate as constant over a timestep uses the condition ΔΛ(t)Δt ≪ 1.This assumption is more lenient than the Λ(t)Δt ≪ 1 condition used by rejection sampling.Source: F. Non-Markovian processes · p. 11
• Updating all time-dependent rates at every timestep is slow, while a first-order cumulant expansion is faster but is valid only when M(t) ≫ 1.The approximation can break down near the beginning and end of an SIR process and near the SIS epidemic threshold.Source: F. Non-Markovian processes · p. 11
• The control parameter ϵ sets rate-refresh timing relative to each process’s mean waiting time, balancing precision against simulation speed.For Weibull recovery times, the mean waiting time is Γ(1 + 1/γ)/µ0.Source: F. Non-Markovian processes · p. 13
• With ϵ = 0.1, the Weibull-recovery SIR example has errors of no more than a few percent with almost optimal computation time.Smaller ϵ improves precision, while larger ϵ increases speed.Source: F. Non-Markovian processes · p. 13

III. DISCUSSION

• For Poisson processes, the temporal Gillespie algorithm is stochastically exact and up to multiple orders of magnitude faster than current algorithms for time-varying networks.Source: §III · p. 13
• For non-Markovian processes, a control parameter lets users choose the desired accuracy and simulation speed.Source: §III · p. 14
• The algorithm is illustrated on compartmental contagion models but may also apply to diffusion-like processes or systems without a network description.Source: §III · p. 14

METHODS

• The supporting information begins with a subsection listing notation used in the article.Source: METHODS · p. 14
• A subsection details how Monte Carlo simulations were performed.Source: METHODS · p. 14
• A subsection provides pseudocode for applying the temporal Gillespie algorithm to specific contagion processes on time-varying networks.Source: METHODS · p. 14
• A final subsection gives pseudocode for optimizing the algorithm on empirical networks by removing obsolete contacts.Source: METHODS · p. 14

Notation

• Table II lists notation pertaining to the temporal Gillespie algorithm.Source: Notation · p. 14
• Table III lists notation for time-varying networks and compartmental contagion processes.Source: Notation · p. 14
• The notation is organized into separate tables for algorithm terms and network and contagion terms.Source: Notation · p. 14

• Speed tests ran sequentially on an Intel Core i7-3687U laptop with 8 GB RAM and a 256 GB SSD.The machine had a dual-core, four-thread CPU and 8 GB of 1,600 MHz DDR3 SDRAM.Source: performed · p. 14
• -O2 was fastest among the tested -O2, -O3, and -Ofast settings for both algorithms.Comparisons were compiled with TDM GCC 64-bit.Source: performed · p. 14
• Each data point used 100–10,000 independent realizations depending on µ∆t; SIR stopped at I = 0, while SIS stopped at 20/(µ∆t) steps or I = 0, whichever came first.Simulations on empirical contact data looped the datasets when necessary.Source: performed · p. 14

• The network is supplied as contact lists indexed by time step, while a state vector records each node as S, I, or R.Time is normalized to ∆t = 1, so code parameters beta and mu correspond to β∆t and µ∆t.Source: Algorithms for simulating specific… · p. 14
• In SIR, infection occurs on contact at rate β and infectious nodes recover spontaneously at rate µ, entering an immune recovered state.The model has susceptible, infectious, and recovered states.Source: Algorithms for simulating specific… · p. 14
• In SIS, infection occurs on contact at rate β and infectious nodes spontaneously return to susceptibility at rate µ.The implementation modifies the SIR pseudocode by removing two lines and replacing the recovery-state update.Source: Algorithms for simulating specific… · p. 14
• C++ implementations are provided for homogeneous and heterogeneous populations.The code covers both SIR and SIS models.Source: Algorithms for simulating specific… · p. 15
• For non-Markovian SIR, recovery times are Weibull distributed and the recovery rate depends on time since infection, while infection remains contact-dependent at constant rate β.The pseudocode treatment also covers recovery times following a general distribution.Source: Algorithms for simulating specific… · pp. 14–15
• The non-constant-rate implementation tracks infection times and the exact time of the last transition; µ(t-t_inf[m]) returns the current recovery rate for node m.µ_avg and epsilon are used in the approximate simulation scheme.Source: Algorithms for simulating specific… · p. 15

empirical networks

• Removing contacts to nodes already in an absorbing state can further speed temporal Gillespie when finite empirical contact data are looped, as in SIR.The procedure replaces the contact-list handling lines of the SIR pseudocode.Source: empirical networks · p. 15

Supplementary Figures

• On high-school face-to-face contact data, the supplementary comparison covers SIR and SIS trajectories and outcome distributions across two rate settings, with 1,000,000 simulations.The outcomes include SIR final epidemic size and SIS stationary infected-node counts.Source: Figure 1 · p. 22
• For high-transition-probability synthetic activity-driven networks, the comparison includes SIR and SIS trajectories and final or stationary-state distributions.Each setting used 1,000,000 simulations on a 100-node network.Source: Figure 2 · p. 23
• In non-Markovian SIR with Weibull recovery times, rejection-sampling outcomes approach temporal Gillespie outcomes as β∆t and µ∆t decrease.The temporal Gillespie precision parameter was set to ϵ = 0 (quasi-exact).Source: Figure 3 · p. 24