Temporal Gillespie Algorithm: Fast Simulation of Contagion Processes on Time-Varying Networks

Stochastic simulations are one of the cornerstones of the analysis of dynamical processes on complex networks, and are often the only accessible way to explore their behavior. The development of fast algorithms is paramount to allow large-scale simulations. The Gillespie algorithm can be used for fast simulation of stochastic processes, and variants of it have been applied to simulate dynamical processes on static networks. However, its adaptation to temporal networks remains non-trivial. We here present a temporal Gillespie algorithm that solves this problem. Our method is applicable to general Poisson (constant-rate) processes on temporal networks, stochastically exact, and up to multiple orders of magnitude faster than traditional simulation schemes based on rejection sampling. We also show how it can be extended to simulate non-Markovian processes. The algorithm is easily applicable in practice, and as an illustration we detail how to simulate both Poissonian and non-Markovian models of epidemic spreading. Namely, we provide pseudocode and its implementation in C++ for simulating the paradigmatic Susceptible-Infected-Susceptible and Susceptible-Infected-Recovered models and a Susceptible-Infected-Recovered model with non-constant recovery rates. For empirical networks, the temporal Gillespie algorithm is here typically from 10 to 100 times faster than rejection sampling.


Author Summary
When studying how e.g. diseases spread in a population, intermittent contacts taking place between individuals-through which the infection spreads-are best described by a timevarying network. This object captures both their complex structure and dynamics, which crucially affect spreading in the population. The dynamical process in question is then usually studied by simulating it on the time-varying network representing the population. Such simulations are usually time-consuming, especially when they require exploration of different parameter values. We here show how to adapt an algorithm originally proposed in 1976 to simulate chemical reactions-the Gillespie algorithm-to speed up such simulations. Instead of checking at each time-step if each possible reaction takes place, as traditional rejection sampling algorithms do, the Gillespie algorithm determines what reaction takes place next

Introduction
Networks have emerged as a natural description of complex systems and their dynamics [1], notably in the case of spreading phenomena, such as social contagion, rumor and information spreading, or epidemics [1][2][3]. The dynamics of contagion processes occurring on a network are usually complex, and analytical results are attainable only in special cases [3,4]. Furthermore, such results almost systematically involve approximations [3,4]. Numerical studies based on stochastic simulations are therefore necessary, both to verify analytical approximations, and to study the majority of cases for which no analytical results exist. The development of fast algorithms is thus important for the characterization of contagion phenomena, and for large-scale applications such as simulations of world-wide epidemics [2,5].
The Doob-Gillespie algorithm [6][7][8][9][10][11] (also known as Gillespie's Stochastic Simulation Algorithm-SSA or Gillespie's direct method), originally proposed by David Kendall in 1950 for simulating birth-death processes and made popular by Daniel Gillespie in 1976 for the simulation of coupled chemical reactions, offers an elegant way to speed up such simulations by doing away with the many rejected trials of traditional Monte Carlo methods. Instead of checking at each time-step if each possible reaction takes place, as rejection sampling algorithms do, the Gillespie algorithm draws directly the time elapsed until the next reaction takes place and what reaction takes place at that time. It is readily adapted to the simulation of Poisson processes on static networks [12][13][14][15][16] and can be generalized to non-Markovian processes [17].
Systems in which spreading processes take place, e.g., social, technological, infrastructural, or ecological systems, are not static though. Individuals create and break contacts at time-scales comparable to the time-scales of such processes [18][19][20], and the dynamics of the networks themselves thus profoundly affect dynamical processes taking place on top of them [21][22][23][24][25][26][27]. This means that one needs to take the network's dynamics into account, e.g., by representing it as a time-varying network (also known as a time-varying graph, temporal network, or dynamical network) [28]. The dynamical nature of time-varying networks makes the adaptation of the Gillespie algorithm to such systems non-trivial.
The main difficulty in adapting the Gillespie algorithm to time-varying networks is taking into account the variation of the set of possible transitions and of their rates at each time step. We show that by normalizing time by the instantaneous cumulative transition rate, we can construct a temporal Gillespie algorithm that is applicable to Poisson (constant rate) processes on time-varying networks. We give pseudocode and C++ implementations for its application to simulate the paradigmatic Susceptible-Infected-Susceptible (SIS) and Susceptible-Infected-Recovered (SIR) models of epidemic spreading, for both homogeneous and heterogeneous [29] populations. We verify the accuracy of the temporal Gillespie algorithm numerically by comparison with a classical rejection sampling algorithm, and we show that it is up to * 500 times faster for the processes and the parameter ranges investigated here.
While Poissonian models are of widespread use, real contagion phenomena show memory effects, i.e., they are non-Markovian. Notably, for realistic infectious diseases, the rate at which an infected individual recovers is not constant over time [30,31]. Social contagion may also show memory effects, e.g., one may be more (or less) prone to adopt an idea the more times one has been exposed to it. To treat this larger class of models, we show how the temporal Gillespie algorithm can be extended to non-Markovian processes. We give in particular an algorithm for simulating SIR epidemic models with non-constant recovery rates.

Results
The following subsections present the main results of the article. Section 1: "Stochastic processes on time-varying networks" defines the stochastic processes which can be simulated using the temporal Gillespie algorithm, and describes the class of compartmental models for contagion phenomena on networks-the class we will use in examples throughout this paper. Section 2: "Rejection sampling for Monte Carlo simulations" gives a quick overview of the traditional rejection sampling algorithms. Section 3: "Gillespie algorithm on static networks" outlines a derivation of the static Gillespie algorithm. Section 4: "Temporal Gillespie algorithm" derives the temporal Gillespie algorithm for Poisson (constant-rate) processes. In Section 5: "Comparison of Gillespie and rejection sampling algorithms" we validate the temporal Gillespie algorithm through numerical comparison with a rejection sampling algorithm; we also compare their speeds for simulating SIR and SIS processes on both synthetic and empirical time-varying networks. Section 6: "Non-Markovian processes" shows how the temporal Gillespie algorithm can be extended to simulate non-Markovian processes; the approach is verified numerically and the speed of the non-Markovian temporal Gillespie algorithm is compared to rejection sampling.
Tables listing the notation used in the manuscript, details on how Monte Carlo simulations were performed, and pseudocode for application of the temporal Gillespie algorithm are given in the Methods section.

Stochastic processes on time-varying networks
We define in this section the type of stochastic processes for which the temporal Gillespie algorithm can be applied. At the time of writing, the main domain of application of the algorithm is the class of compartmental models for contagion processes on time-varying networks, which we introduce below. For definiteness, algorithms detailing the application of the temporal Gillespie algorithm will concern this class of stochastic processes.
In general, we consider a system whose dynamics can be described by a set of stochastic transition events. We assume that the system can be modeled at any point in time by a set, O(t), of M(t) independent stochastic processes m, which we term transition processes; the rate at which the transition m takes place is denoted λ m . The set O(t) thus defines the possible transition events at time t and in general changes over time, depending on both external factors and the evolution of the system itself; the number of possible transitions, M(t), thus also generally changes over time, while λ m may or may not vary over time. For the classic "static" Gillespie algorithm to be applicable, O(t) is allowed to change only when a transition (or chemical reaction in the context of Gillespie's original article) takes place. For processes taking place on time-varying networks, the medium of the process-the network-also changes with time. As these changes may allow or forbid transitions, O(t) is not only modified by every reaction, but also by every change in the network. This is the domain of the temporal Gillespie algorithm, which only requires that the number of points in which O(t) changes be finite over a finite time-interval [32].
The assumption that the transition processes are independent is essential to the validity of the Gillespie algorithm, as it allows the calculation of the distribution of waiting times between consecutive transitions. This assumption is not overly restrictive, as it only requires a transition process to be independent of the evolution of the other simultaneous transition processes. A transition process may depend on all earlier transitions, and the current and past states of all nodes. As such, Gillespie algorithms can notably be applied to models of cooperative infections and other non-linear processes such as threshold models [17], and has even been applied to model the dynamics of ant battles [33].
Compartmental models of contagion. In a network-based description of the population in which a contagion process takes place, an individual is modeled as a node i (Fig 1A). A contact between two individuals taking place at time t is represented by an edge (i, j) t in the graph describing the population at the instant t (Fig 1A). In a compartmental model, each node is in a certain state, which belongs to a fixed, finite set of q different states (compartments) [3]. A random variable x i ðtÞ 2 fX 1 ; X 2 ; . . . ; X q g specifies the state of the node i at time t (i.e. to which compartment it belongs). Nodes may stochastically transition between states, governed by the set O(t) of transition processes. One is usually interested in the evolution of the number of nodes in each state, which we denote X 1 , X 2 , . . ., X q .
As an example, we consider the classic SIR model of epidemic spreading with constant transition rates in a homogeneous population (rates are the same for all individuals) ( Fig 1B). Here nodes can be in one of three states: susceptible, infected, and recovered, fS; I ; Rg. Two different transition types let nodes change state: (i) a node i in the S state switches to the I state with rate state with rate k I ðtÞb (an S ! I reaction), where k I ðtÞ is the number of contacts i has with nodes in the I state at time t ( Fig 1A); (ii) a node i in the I state switches to the R state at rate μ (an I ! R reaction).
In general, the transition processes can be divided into three types: 1. spontaneous transitions, which only depend on the current state of the node, x i (t) (Fig 1C) -e.g. an infected node recovers spontaneously in the SIR model ( Fig 1B); 2. contact-dependent transitions, which may take place only when the node i is in contact with other nodes in a given state; it thus depends on the states x j of the nodes j currently in contact with i ( Fig 1D)-e.g. a susceptible node may be infected upon contact with an infectious node in the SIR model ( Fig 1B).
3. mixed transitions, which take place spontaneously, but may depend on the node's past and current contacts (Fig 1E)-e.g. in rumor spreading, an individual may learn on his own that a rumor is false (spontaneous) or may be convinced by another individual who knows the rumor is false (contact-dependent).
This division is important for practical application of the temporal Gillespie algorithm as transition processes of type a need only be updated after a transition has taken place, and processes of type c need only be updated whenever a relevant contact takes place, but not at each time-step. Using this distinction is crucial to obtaining the large speed-increase that the temporal Gillespie algorithm offers over rejection sampling, as discussed below (Secs. 4: "Temporal Gillespie algorithm" and 5: "Comparison of Gillespie and rejection sampling algorithms").

Rejection sampling for Monte Carlo simulations
A straightforward way to simulate a stochastic process is to use a rejection sampling algorithm, akin to the classical Metropolis algorithm. Here one divides the time-axis in small time-steps Δt, where Δt should be chosen sufficiently small such that this discretization does not influence the outcome of the process significantly; on time-varying networks, the choice of Δt often comes naturally as the time-resolution at which the network is measured or simulated (Fig 1A).
At each time-step t = 0, Δt, 2Δt, . . ., we test whether each possible transition m 2 O(t) takes place or not. In practice this is done by drawing a random number r m that is uniformly distributed on [0, 1) for each m and comparing it to λ m Δt: if r m < λ m Δt the reaction takes place, if r m ! λ m Δt nothing happens [Fig 2 (Transitions)]. (Note that one should technically compare r m to 1 − exp(λ m Δt) to ensure that λ m defines a proper transition rate for finite Δt. However, the two procedures are equivalent in the limit Δt ! 0.) From the design of the rejection sampling algorithm we see that the proportion of trials that are rejected is equal to a weighted average over {1 − λ m Δt} m . Thus, since we require λ m Δt ( 1 Nodes' colors correspond to their current state; edges denote current contacts between nodes; edge colors correspond to: black: no contagion may take place over the edge, red: contagion takes place during the present time-step, and red-to-blue gradient: contagion is possible but does not take place. (B) Example: reaction types in the SIR model. (C) Spontaneous reaction: a node i may spontaneously transition from its current state x i to x 0 i with rate λ m . (D) Contact-dependent reaction: a node i may transition from its current state x i to x 0 i with rate λ m upon contact with the node j in state x j . (E) Mixed transition: a node i may spontaneously transition from its current state x i to another state, x 0 i with rate λ m ; contact with another node j, in state x j , may alter the transition rate of m, l m ! l 0 m . After the contact (i, j) t ends, the transition rate may revert to λ m , remain unchanged, or change to a third value.  . From the state of the network at the first time-step, the set of possible transitions Ω(0) is found (Transitions), and from this the cumulative transition rate Λ(0) is calculated. The integrated cumulative transition rate, LðDt; 0Þ ¼ Lð0ÞDt is compared to t 0 1 . If, as in the present example, Lð0ÞDt < t 0 1 the algorithm is advanced to the next time-step. (B) and (C) The set of possible transitions Ω(t n ) and the cumulative transition rate Λ(t n ) is updated at each following time-step n; if Lðt n ; 0Þ ¼ Dt P nÀ1 l¼0 Lðt l Þ is still smaller than t 0 1 , the algorithm is advanced to the next time-step. (D) During the first time-step n** where Lðt n ÃÃ ; 0Þ ! t 0 1 , a in order to avoid discretization errors, the vast majority of trials are rejected and the rejection sampling algorithm is computationally inefficient.

Gillespie algorithm on static networks
The Gillespie algorithm lets us perform stochastically exact Monte Carlo simulations without having to reject trials. For Poisson processes on static networks, it works by recognizing that the waiting time between two consecutive transitions is exponentially distributed, and that each transition happens with a probability that is proportional to its rate. Specifically, the (survival) probability that the transition m has not taken place after a time τ since the last transition event is Since each transition takes place independently, the probability that no event takes place during the interval τ since the last event is where L ¼ P M m¼1 l m is the cumulative transition rate. The above result is obtained by using the fact that while O and M do depend on t, they only change when an event takes place and not in-between. They can thus be treated as constant for the purpose of calculating the waiting time between events. The distribution of the waiting times τ is then given by the probability density p(τ) = Λe − Λτ , while the probability density for the reaction m being the next reaction that takes place and that it takes place after exactly time τ is equal to p m (τ) = λ m e − Λτ The static Gillespie algorithm thus consists in drawing the waiting time τ* Exp (Λ) until the next transition and then drawing which transition m takes place with probability π m = λ m /Λ. [Here τ* Exp (Λ) is short for: τ is exponentially distributed with rate Λ.]

Temporal Gillespie algorithm
For processes taking place on time-varying networks however, the set of transition process, O (t), changes with time independently of the transition events, e.g., for the case of an SIR process nodes may become infected only when in contact with an infected individual ( Fig 1A). This means that the survival probability does not reduce to a simple exponential as in Eq (1); it is instead given by where t Ã is the time at which the last transition took place, t ÃÃ = t Ã + τ is the time when the next transition takes place, and I m (t) is an indicator function that is equal to one when the process m may take place, e.g., when two given nodes are in contact, and zero when m may not take place. The meaning of I m is exemplified in Fig 1A: the node i may be infected by the infectious node j only when the two nodes are in contact; if we let m denote this transition process, I m (t) is then one for t = Δt, 3Δt, 4Δt and zero for t = 0, 2Δt.
transition takes place. The exact time of this transition, t**, is given by Eq (12) and the transition m that takes place is chosen among the possible transitions in the given time-step with probability proportional to its transition rate λ m [Transitions and Eq (10)]. (E) The transition changes the system (Network) and consequently the set of possible transitions, Ω(t**), (Transitions); thus Ω(t**) and Λ(t**) must be updated, which in turn changes the remainder of Lðt n ÃÃ ; 0Þ (Normalized time). A new normalized waiting time is then drawn, t 0 2 $ Expð1Þ; if Lðt n ÃÃ ; 0Þ < t 0 1 þ t 0 2 , no further transitions takes place during the time-step and the algorithm is advanced to the next time-step (note that multiple transitions may occur during the same time-step). (F) The above procedure is reiterated. Note that for processes taking place on adaptive time-varying networks, whose changes only depend on the process itself, I m (t) only changes when a transition takes place and Eq (3) reduces to Eq (1). This means that from the point of view of the algorithm, such networks are effectively static and the classic "static" Gillespie algorithm may simply be used there [14,16].
We now consider the general case where O(t) may change independently of the processes evolving on the network (described in Sec. 1: "Stochastic processes on time-varying networks"). Using, as in the previous section, that transition processes are independent, we can write the probability that no event takes place during an interval τ (the waiting time survival function): where O denotes the set of all possible transitions (transition processes) on the interval between two transition events, (t Ã , t ÃÃ ], i.e., O is the union over O(t) for t 2 (t Ã , t ÃÃ ], and M is the total number of transition processes on the same interval (the size of O). We switch the sum and the integral in Eq (4) to obtain Finally, using that I m (t) = 0 for all m = 2 O(t), we may write where is the cumulative transition rate at time t.
The dynamics of empirical time-varying networks is highly intermittent and we cannot describe O(t) analytically. This means that we cannot perform the integral of Eq (6) to find the waiting time distribution directly. We may instead normalize time by the instantaneous cumulative transition rate, Λ(t): We define a unitless normalized waiting time between two consecutive transitions, τ 0 , as i.e., equal to the cumulative transition rate integrated over (t Ã , t ÃÃ ]. The survival function of τ 0 has the following simple form: The time t ÃÃ when a new transition takes place is given implicitly by Lðt ÃÃ ; t Ã Þ ¼ t 0 , while the probability that m is the transition that takes place at time t = t ÃÃ is given by: This lets us define a Gillespie-type algorithm for time-varying networks by first drawing a normalized waiting time τ 0 until the next event from a standard exponential distribution [i.e. with unit rate, τ 0 * Exp (1)], and second, solving Lðt; t Ã Þ ¼ t 0 numerically to find t ÃÃ . In practice, since Λ(t) only changes when a transition takes place or at t n = nΔt with n 2 N, we need only compare τ 0 to for each time-step n (Fig 2A-2C). Here n Ã is the time-step during which the last transition took place, and Λ(t Ã ) is the cumulative transition rate at t Ã , immediately after the last transition has taken place. The first term of Eq (11) is the cumulative transition rate integrated over the remainder of the n Ã th time-step left after the last transition; the second term is equal to Lðt nþ1 ; t n Ã þ1 Þ. A new transition takes place during the time-step n ÃÃ where Lðt n ÃÃ þ1 ; t Ã Þ ! t 0 ( Fig  2D); the precise time of this new transition is the reaction m that takes place is drawn with probability given by Eq (10) (Fig 2D). We then update O and Λ to O(t ÃÃ ) and Λ(t ÃÃ ) (Fig 2E), draw a new waiting time, τ 0 * Exp (1), and reiterate the above procedure ( Fig 2F). The algorithm can be implemented for contagion processes on time-varying networks as follows (see Methods for pseudocode for specific contagion models and S1 Files for implementation in C++): 1. Draw a normalized waiting time until the first event from a standard exponential distribution, τ 0 * Exp (1) (Fig 2A).
2. At each time-step t n = nΔt, with n = 0, 1, 2, . . ., let O O(t n ) and Λ Λ(t n ); here, only contact-dependent processes (type b, Sec. 1: "Stochastic processes on time-varying networks") and mixed (type c) processes that depend on contacts taking place at t n or t n − 1 need to be updated-an important point, as it lets the temporal Gillespie algorithm be much faster than rejection sampling (see Discussion in the following section). Then, compare τ 0 to ΛΔt: if ΛΔt τ 0 Subtract ΛΔt from τ 0 , continue to next time-step and repeat 2 (Fig 2A-2C) [34].
if ΛΔt > τ 0 Let the reaction m take place, chosen from O with probability π m = λ m /Λ. The fraction that is left of the time-step when the transition takes place is ξ = 1 − τ 0 /(ΛΔt) and the precise time of the transition is t ÃÃ = t n + τ 0 /Λ (Fig 2D and 2F). Next, update O and Λ ( Fig  2E); this time all transition processes should be updated, as spontaneous processes (type a) may change, emerge, or disappear when a transition takes place. Then: (a) draw a new normalized waiting time, τ 0 * Exp (1) ( Fig 2F); if τ 0 ! ξΛΔt subtract ξΛΔt from τ 0 , continue to the next time-step and repeat 2 (Fig 2F).
if τ 0 < ξΛΔt Another transition takes place during the present time-step (at time t ÃÃÃ = t ÃÃ + τ 0 /Λ, where t ÃÃ is the time of the last transition during the same time-step): choose m from O with probability π m = λ m /Λ; let ξ ! ξ − τ 0 /ΛΔt, and update O and Λ. Repeat a) and b).
By construction, the above procedure produces realizations of a stochastic process for which the waiting times for each transition follow exactly their correct distributions. The temporal Gillespie algorithm is thus what we term stochastically exact: all distributions and moments of a stochastic process evolving on a time-varying network obtained through Monte Carlo simulations converge to their exact values. Rejection based sampling algorithms are stochastically exact only in the limit λ m Δt ! 0.
A large literature exists on the related problem of simulating coupled chemical reactions under externally changing conditions (e.g., time-varying temperature or volume) [35][36][37][38][39][40]. Most of these methods consider only external perturbations that can be described by an analytical expression. In this case the problem reduces to that of defining a static, yet non-Markovian, algorithm. Some methods, and notably the modified next reaction method developed by Anderson [37], can be adapted to a completely general form of the external driving and thus, in principle, to simulate dynamical processes taking place on time-varying networks. These methods are based on a scheme that is conceptually similar to Gillespie's direct algorithm, the next reaction method, proposed by Gibson and Bruck [35]. The next reaction method draws a waiting time for each reaction individually and chooses the next reaction that happens as the one with the shortest corresponding waiting time. It then updates the remaining waiting times, draws new waiting times (if applicable), and reiterates. To generalize the next reaction method to processes with non-exponential waiting times, Anderson introduced the concept of the internal time for each transition process [37]. In the notation used in the present article it is defined as T m ðtÞ ¼ R t 0 I m ðtÞl m dt and is thus equivalent to the normalized time, Lðt; 0Þ, only for an individual transition process.
By construction, the next reaction method needs to draw only one random number per transition event, where the Gillespie algorithm draws two. However, this reduction in the number of required random variables comes at a price: one must draw a random number for each individual transition process and keep track of, compare, and update each of the individual waiting times. For chemical reactions, where the number of different chemical reactions is small (it scales with the number of chemical species), this tradeoff favors the next reaction method. However, for contagion processes on networks, each individual is unique (if not intrinsically, at least due to its position in the network). The number transition processes thus scales with the number of nodes and contacts, which favors the Gillespie algorithm as it does not need to keep track of each of them individually [17].
On time-varying networks (or for time-varying external driving) one must furthermore update relevant internal times each time the network structure (external conditions) changes in the next reaction method. Chemically reacting systems are usually close to being adiabatic, i.e., the external driving changes slowly compared to the time-scales of chemical reactions. Thus, the additional overhead related to updating individual internal times is practically negligible. However, the dynamics of temporal networks is highly intermittent and the time-scale of network change is typically smaller than the time-scales of relevant dynamical processes. Here one must thus update the internal times many times between each transition event, inducing a substantial overhead. Since the temporal Gillespie algorithm operates with a single global normalized waiting time, it handles these updates more efficiently.
Finally, the modified next reaction method may in principle be extended to non-Markovian processes taking place on time-varying networks (as treated in Sec. 6: "Non-Markovian processes" using the temporal Gillespie algorithm). However, such an approach would, for each single transition, require solving numerically Eq. (13) of [37] for the internal waiting time of each individual transition process, taking into account the time-varying network structure, finding the shortest corresponding waiting time in real time, and then updating the internal waiting times of all the other reactions, rendering the next reaction method even more inefficient in this general case.

Comparison of Gillespie and rejection sampling algorithms
Numerical validation. We compare the outcome of SIR and SIS processes on activitydriven time-varying networks [41] simulated using the temporal Gillespie algorithm to the outcome of simulations using traditional rejection sampling. For sufficiently small λ m Δt, the outcomes are indistinguishable (Fig 3, see also S1 Fig for an empirical network of face-to-face contacts in a high school), confirming the validity of the temporal Gillespie algorithm. Note that rejection sampling is only expected to be accurate for λ m Δt ( 1, while the temporal Gillespie algorithm is stochastically exact for all λ m Δt; the results of the two algorithms thus differ when the assumption λ m Δt ( 1 does not hold (S2 Fig). Comparison of simulation speed. Next, we compare the speeds of the temporal Gillespie and the rejection sampling algorithms for SIR and SIS processes (see Methods for details on how simulations were performed).   Table 1). The speed gain is higher for larger systems (compare N = 1 000 to N = 100 in Fig 4) We also see that the speed gain is larger the sparser the network is. This is because the calculation of the contacts between susceptible and infected nodes at each time-step, necessary to determine the possible S ! I transitions, is the performance limiting step of the temporal Gillespie algorithm (see below). In the extreme case of a contagion model where all transitions are contact-dependent (type b, Sec. 1: "Stochastic processes on time-varying networks"), such as the classic Maki-Thompson model of rumor spreading [42], the temporal Gillespie algorithm is approximately a factor two faster than the rejection sampling algorithm.
Expected time complexity of the algorithms. We may gain insight into the performance of the algorithms by considering their time-complexity, i.e., how their running time scales with the input parameters of the simulated system. Since the algorithms are used for Monte Carlo simulations, it is most interesting to consider the expected complexity given a set of parameters, i.e., the mean running time of an algorithm averaged over an ensemble of simulations, not the worst-case complexity which is usually considered for deterministic algorithms.
The expected running time of the rejection sampling algorithm scales as where OðxÞ denotes a term that is of order x, EðtÞ ¼ NkðtÞ=2 is the mean number of contacts  Gillespie algorithm is given by where QðtÞis the mean number of transitions that take place per time-step. The first term of the r.h.s. of Eqs (13) and (14) correspond to the time needed for looking through the set of contacts at each time-step to determine the set of possible infections and are thus similar for the rejection sampling and temporal Gillespie algorithms (with the temporal Gillespie algorithm incurring a small additional overhead related to calculating the cumulative transition rate and keeping track of of the normalized waiting time left till the next transition). For rejection sampling [Eq (13)], the second term corresponds to the determination of whether each of the possible transitions takes place at each time-step; for the temporal Gillespie algorithm [Eq (14)], the second term corresponds to drawing inter-event waiting times and which transitions that take place. For the SIR and SIS processes considered above,   (14)]; this explains why the difference in performance decreases with the mean instantaneous degree of the network. This also hints at how we may improve the speed of the temporal Gillespie algorithm: by rendering the identification of relevant contacts during each time-step faster. One such approach which may be applied to processes with an absorbing state (e.g. an R state) is explored below. Improving performance by removing obsolete contacts. Empirical networks describing human contact differ from simulated networks in a number of ways. For example, their structure and dynamics are more complex [25,[46][47][48][49] but perhaps most importantly in the perspective of optimizing simulations, they are of finite length. One is often interested in long-time behavior or slowly evolving processes compared to the length of available data. To overcome this limitation, one usually loops over the data set. This means that if a node enters an inactive absorbing state such as the recovered (R) state in the SIR model, one may remove all following contacts to this node from the data, thus reducing the number of contacts that one must go through during the following loop. Furthermore, since the I ! R transition is independent of the network, one may also remove all contacts between two infected nodes.
Pseudocode for an algorithm that removes obsolete contacts is given in Methods and C++ code is given in S1 Files. Fig 5A compares the speed gain of the temporal Gillespie algorithm relative to rejection sampling with and without contact removal for simulations of a constant-rate SIR process on empirical networks of face-to-face contacts (Table 1). Depending on the parameters of the simulated process, removing obsolete contacts may induce both a significant gain or loss in speed; for processes that are fast compared to the length of the data set, the data is not repeated or only repeated few times during a simulation and the additional overhead involved in identifying and removing the obsolete contacts renders the algorithm slower; for slow processes the data is looped many times and removing the obsolete contacts makes the algorithm faster. Fig 5A suggests an empirically determined rule-ofthumb: if the slowest time-scale of the simulated process (here * 1/μ) is longer than the length of the data, T, removing obsolete contacts pays off, if it is shorter, one should not remove obsolete contacts.
Slow network dynamics. For time-varying networks of face-to-face contacts, which are relevant for simulating epidemic spreading in a population, network dynamics are typically much faster than the time-scales of the dynamical process that is simulated (compare the 20 s timeresolution of the empirical data of Table 1 to typical 1/β * 1 hour and 1/μ * 1 week for flu-like diseases). In the opposite case, i.e., if the network evolves much slower than the dynamical process, the temporal Gillespie algorithm simply works like a static Gillespie algorithm in-between changes in the network structure while taking the changes changes into account exactly when they occur. The performance of the temporal Gillespie algorithm then approaches that of a static Gillespie algorithm in this case. Note that since QðtÞ ) EðtÞin this limit, the second term dominates in Eq (14), which means that the speed of the algorithm is limited by the selection of waiting times and transitions that take place, and care should be taken to optimize these steps, e.g., by organizing the transition processes in a heap or a priority queue [37]. Note finally that to obtain reliable results using a rejection sampling algorithm one must use a time-step for simulations Δt RS which is much smaller than the time-step Δt of network change. Thus the expected time complexity of rejection sampling scales with Δt/Δt RS n simu ) n simu in this case.

Non-Markovian processes
For real-world contagion processes, transition rates are typically not constant but in general depend on the history of the process [30,31]. Such processes are termed non-Markovian. The survival probability for a single non-Markovian transition process taking place on a time-varying network is given by: Here F ðmÞ t is a filtration for the process m, i.e., all information relevant to the transition process available up to and including time t; typically, F ðmÞ t will be its starting time and relevant contacts that have taken place since. As above, t Ã is the time of the last transition and t ÃÃ = t Ã + τ is the time of the next. [Note that since λ m now depends explicitly on t, we may absorb I m in λ m ; however, to underscore the analogy with the Poissonian case, we keep the factor I m explicitly in Eq (15). ] We use again that the transition processes are independent, to write the waiting time survival probability: and where F t is the union over F ðmÞ t for m 2 O. For a static network, Eq (6) reduces to the result found in [17]. This can be seen by noting that M(t) = M and O(t) = O are then constant, and thus that l m ðt; F ðmÞ t Þ ¼ À½dS m ðt; F ðmÞ t Þ=dt=S m ðt; F ðmÞ t Þ ¼ df ln ½1=S m ðt; F ðmÞ t Þg=dt and S m ðt; F ðmÞ t Þ ¼ S m ðt þ t m ; F ðmÞ t Þ=S m ðt m ; F ðmÞ t Þ, yielding directly Eq. (7) of [17]. As in the Poissonian case (Sec. 4: "Temporal Gillespie algorithm") we define the normalized waiting time, τ 0 , as This gives us the same simple form as above for the survival function of the normalized waiting time, τ 0 , and the probability that m is the transition that takes place at t = t ÃÃ , Until now our approach and results are entirely equivalent to the Poissonian case considered above. However, since λ m (t) in general depend continuously on time, the transition time t ÃÃ is not simply found by linear interpolation as in Eq (12). Instead, one would need to solve the implicit equation Lðt ÃÃ ; t Ã Þ ¼ t 0 numerically to find t ÃÃ exactly. To keep things simple and speed up calculations, we may approximate Λ(t) as constant over a time-step. This assumes that ΔΛ(t)Δt ( 1, where ΔΛ(t) is the change of Λ(t) during a single time-step. It is a more lenient assumption than the assumption that Λ(t)Δt ( 1 which rejection sampling relies on, as can be seen by noting that in general ΔΛ(t)/Λ(t) ( 1. The same assumption also lets us calculate Lðt nþ1 ; t Ã Þ as in the Poissonian case: and the time, t ÃÃ , at which the next transition takes place: Using the above equations, we can now construct a temporal Gillespie algorithm for non-Markovian processes. This algorithm updates all λ m (t) that depend on time at each time-step, where for the Poissonian case we only had to initialize new processes, i.e., contact-dependent processes (type b and c, Sec. 1: "Stochastic processes on time-varying networks"). This means the algorithm is only roughly a factor two faster than rejection sampling [compare dotted lines ( = 0) in Fig 6]. To speed up the algorithm, we may employ a first-order cumulant expansion of Sðt; F t Þ around τ = 0, as proposed in [17,38] for static non-Markovian Gillespie algorithms. It consists in approximating Temporal Gillespie Algorithm l m ðt; F ðmÞ t Þ by the constant l m ðt Ã ; F ðmÞ t Þ for t Ã < t < t ÃÃ and gives a considerable speed increase of the algorithm [full ( ! 1) in Fig 6]. However, the approximation is only valid when M(t) ) 1 [43], which is not always the case for contagion processes. Notably, at the beginning and end of an SIR process, and near the epidemic threshold for an SIS process, M is small and the approximation breaks down; the approximate algorithm for example overestimates the peak number of infected nodes in a SIR process with recovery rates that increase over time [compare full black line ( ! 1) to the quasi-exact full red line ( = 0) in Fig 7A]. An intermediate approach, which works when the number of transition processes is small, but is not too slow to be of practical relevance, is needed. We propose one such approach below [44].
Efficient non-Markovian temporal Gillespie algorithm. As discussed above, we neither want to update all transition rates at each time-step as this makes the temporal Gillespie algorithm slow, nor do we want to only update them when a transition event takes place as this makes the algorithm inaccurate.
An intermediate approach is found by looking at the relevant physical time-scales of the transition processes: the average waiting time before they take place, hτ (m) i. If the time elapsed since we last updated λ m (t) is small compared to hτ (m) i, we do not make a large error by treating it as constant over the interval; however, if the elapsed time is comparable to or larger than hτ (m) i, the error may be considerable. Thus, instead of updating λ m at each time-step, we may update it only after a time t > hτ (m) i has elapsed since it was last updated. Here controls the precision of the algorithm.
Below, we use this approach to simulate a non-Markovian SIR process, where the recovery times of infected nodes follow a Weibull distribution (see Methods for an algorithm written in pseudocode and S1 Files for implementation in C++). The recovery rate of an infected node is here given by where μ 0 sets the scale, t (m) is time when the node was infected, and γ is a shape parameter of the distribution. For γ = 1, we recover the constant-rate Poissonian case with μ(t;t (m) ) = μ 0 . For realistic modeling of infections, γ > 1; here μ(t;t (m) ) is zero at t = t (m) and grows with time. In this case, we thus update the recovery rates μ(t;t (m) ) whenever the time elapsed since a transition last took place exceeds hτ (m) i = Γ(1 + 1/γ)/μ 0 .
The parameter lets us control the precision of the non-Markovian temporal Gillespie algorithm: the smaller is, the more precise the algorithm is, on the other hand, the larger is, the faster the algorithm is (Fig 8). At = 0, the temporal Gillespie algorithm is maximally accurate, but also slowest, corresponding to the quasi-exact approximation that Lðt; F t Þ stays constant over a single time-step. Letting ! 1 corresponds to the first order cumulant expansion of [17], and is the fastest, but least accurate. Intermediate gives intermediate accuracy and speed, and permits one to obtain the desired accuracy without sacrificing performance. In the case of the SIR process with Weibull-distributed recovery times, = 0.1 gives an error of no more than a few percent (Figs 8A-8D and 7)-which is usually acceptable as the discrepancy between model and reality can be expected to be larger-with an almost optimal computation time (Figs 8E and 6).

Discussion
We have presented a fast temporal Gillespie algorithm for simulating stochastic processes on time-varying networks. The temporal Gillespie algorithm is up to multiple orders of magnitude faster than current algorithms for simulating stochastic processes on time-varying networks. For Poisson (constant-rate) processes, where it is stochastically exact, its application is particularly simple. The algorithm is also applicable to non-Markovian processes, where a control parameter lets one choose the desired accuracy and performance in terms of simulation speed. We have shown how to apply it to compartmental models of contagion in human contact networks. The scope of the temporal Gillespie algorithm is more general than this, however, and it may be applied e.g. to diffusion-like processes or systems for which a network description is not appropriate.

Methods
The following four subsections contain supporting information to the manuscript: the first subsection lists notation used in the article (Notation); the second details how Monte Carlo simulations were performed (Details on how Monte Carlo simulations were performed) the third gives pseudocode for application of the temporal Gillespie algorithm to specific contagion processes on time-varying networks (Algorithms for simulating specific contagion models). Finally, in the fourth subsection we give pseudocode for further optimization of the algorithm for empirical networks by removal of obsolete contacts (Removing obsolete contacts for an SIR process on empirical networks).

Notation
Tables 2 and 3 list the notation used in the manuscript. Table 2 gives notation pertaining to the temporal Gillespie algorithm, and Table 3 lists notation pertaining to time-varying networks and compartmental contagion processes.

Details on how Monte Carlo simulations were performed
All simulations for comparing the speed of algorithms were performed sequentially on a HP EliteBook Folio 9470m with a dual-core (4 threads) Intel Core i7-3687U CPU @ 2.10 GHz. The system had 8 GB 1 600 MHz DDR3 SDRAM and a 256 GB mSATA Solid State Drive. Code was compiled with TDM GCC 64 bit using g++ with the optimization option -O2. Speedtests were also performed using -O3 and -Ofast, but -O2 gave the fastest results, both for rejection sampling and temporal Gillespie algorithms.
For SIR processes simulations were run until I = 0; for SIS processes simulations were run for 20/(μΔt) time-steps (as in Fig 3) or until I = 0, whichever happened first. Between 100 and 10 000 independent realizations were performed for each data point depending on μΔt (100 for low μΔt and 10 000 for high μΔt). For simulations on empirical contact data, data sets were looped if necessary.

Algorithms for simulating specific contagion models
We here give pseudocode for the application of the temporal Gillespie algorithm to three specific models: the first subsection treats the Poissonian SIR process, the second treats the Poissonian SIS process, and the third treats a non-Markovian SIR process with recovery times following a general distribution.
We assume in the following that the time-varying network is represented by a list of lists of individual contacts taking place during each time-step. An individual contact, termed Table 2. Notation pertaining to the temporal Gillespie algorithm. The row "First appearance(s)" points to where where the notation is introduced in the Results section.

Symbol
Description First appearance (s) As one may always normalize time by the duration of a time-step, Δt, we have in the following set Δt = 1. Note that beta and mu in the code then corresponds to βΔt and μΔt, respectively. SIR process. The classical SIR model with constant infection and recovery rates is the simplest epidemic model where individuals can gain immunity. As discussed in the main text, nodes may be in one of three states: susceptible (S), infectious (I ), or recovered (R). Two different transition types let the nodes switch state: a spontaneous I ! R transition which takes place with rate μ, and a contact-dependent S ! I transition which takes place with rate β upon contact with an infectious node. Pseudocode 1 shows how the temporal Gillespie algorithm may be implemented for an SIR process on a time-varying contact network. Table 3. Notation pertaining to compartmental contagion models and time-varying networks. The row "First appearance(s)" points to where where the notation is introduced in the Results section. SIS process. In the SIS model, nodes can be in one of two states: susceptible (S) or infectious (I ). As for the SIR model, two different transition types let the nodes switch state: a spontaneous I ! S transition which takes place with rate μ, and a contact-dependent S ! I transition which takes place with rate β upon contact with an infectious node.
The SIS model is implemented in a manner very similar to the SIR model; an implementation can be found by using the code of Pseudocode 1 with lines 07 and 40 removed and line 37 replaced by x[ m] = S. C++ code is given in S1 Files for both homogeneous and heterogeneous populations.
Non-Markovian SIR process. We consider in the main text (Sec. 6: "Non-Markovian processes") an SIR model with non-constant recovery rates; instead of being exponentially distributed (as in the constant-rate SIR model), the individual recovery times, τ (m) , are here Weibull distributed, As for the classical SIR model, nodes may be in one of three states: susceptible (S), infectious (I ), or recovered (R). Two different transition types let the nodes switch state: a contactdependent S ! I transition with constant rate β upon contact with an infectious node, and a spontaneous I ! R transition which takes place with rate μ(t;t (m) ), given by Eq (23). The implementation of the SIR model with non-exponentially distributed waiting times requires some extension of the code for the constant-rate SIR model to account for the heterogeneous and time-dependent recovery rates. To this end, we introduce the following variables: t_inf lists the times at which each node was infected (if applicable); t Ã is the exact time at which the last transition took place; mu is a function of time that is called as mu(t-t_inf[ m] ) to return the instantaneous recovery rate of m at time t; mu_avg is the expected recovery time of an infected node and is used together with the precision control parameter epsilon in the approximate simulation scheme discussed in Section 6: "Non-Markovian processes". Pseudocode 2 shows pseudocode for an implementation of such a SIR model with non-constant recovery rates.

Removing obsolete contacts for an SIR process on empirical networks
When simulations are carried out on data which are looped due to their finite length, the speed of the temporal Gillespie algorithm may be further increased for processes with an absorbing state, such as the SIR process, by removing obsolete contacts to nodes that have entered such a state. Pseudocode 3 shows pseudocode for removing obsolete contacts; its replaces lines 11 to 19 of Pseudocode 1.  The outcome of the rejection sampling algorithm approaches that of the temporal Gillespie algorithm as βΔt and μΔt become smaller. All simulations were performed 100 000 times with the root node chosen at random on an activity driven network consisting of N = 100 nodes, with activities a i = ηz i , where η = 0.1 and z i $ z À3:2 i for z i 2 [0.03,1), and a node formed two contacts when active. Nodes' recovery times followed Eq (20) with γ = 1.5 and the length of a time-step was Δt = 1 s. (PDF) S1 Files. C++ code for implementations of the temporal Gillespie algorithm to examples of epidemic processes on time-varying networks. Specifically, we provide the following programs: (SIR-Poisson-homogeneous.cpp) constant-rate SIR process in a homogeneous population, i.e., same infection/recovery rates for all nodes; (SIR-Poisson-homogeneous-contactRemoval.cpp) constant-rate SIR process in a homogeneous population where obsolete contacts are removed from the contact data as they occur; (SIR-Poisson-heterogeneous.cpp) constant-rate SIR process in a heterogeneous population, i.e., infection/recovery rates may differ between nodes; (SIR-nonMarkovian.cpp) non-Markovian SIR process with Weibull distributed recovery times of individual nodes; (SIS-Poisson-homogeneous.cpp) constant-rate SIS process in a homogeneous population, i.e., same infection/recovery rates for all nodes; (SIS--Poisson-heterogeneous.cpp) constant-rate SIS process in a heterogeneous population, i.e., infection/recovery rates may differ between nodes. (ZIP)