2.1 Jump Processes

Before getting into the nitty-gritty of the Gillespie algorithms, we first explore which types of systems they can be used to simulate. First of all, with the Gillespie algorithms, we are interested in simulating a dynamic system. This can be, for example, epidemic dynamics in a population in which the number of infectious individuals varies over time, or the evolution of the number of crimes in a city, which also varies over time in general. Second, the Gillespie algorithms rely on a predefined and parametrized mathematical model for the system to simulate. Therefore, we must have the set of rules for how the system or the individuals in it change their states. Third, Gillespie algorithms simulate stochastic processes, not deterministic systems. In other words, every time one runs the same model starting from the same initial conditions, the results will generally differ. In contrast, in a deterministic dynamical system, if we specify the model and the initial conditions, the behavior of the model will always be the same. Fourth and last, the Gillespie algorithms simulate processes in which changes in the system are primarily driven by discrete events taking place in continuous time. For example, when a chemical reaction obeying the chemical equation A $+$ B $\to$ C $+$ D happens, one unit each of A and of B are consumed, and one unit each of C and of D are produced. This event is discrete in that we can count the event and say when the event happened, but it can happen at any point in time (i.e., time is not discretized but continuous).

We refer to the class of mathematical models that satisfy these conditions and may be simulated by a Gillespie algorithm as jump processes. In the remainder of this section, we explore these processes more extensively through motivating examples. Then, we introduce some fundamental mathematical definitions and results that the Gillespie algorithms rely on.

2.2 Representing a Population as a Network

Networks are an extensively used abstraction for representing a structured population, and Gillespie algorithms lend themselves naturally to simulate stochastic dynamical processes taking place in networks. In a network representation, each individual in the population corresponds to a node in the network, and edges are drawn between pairs of individuals that directly interact. What constitutes an interaction generally depends on the context. In particular, for the simulation of dynamic processes in the population, the interaction depends on the nature of the process we wish to simulate. For simulating the spread of an infectious disease, for example, a typical type of relevant interaction is physical proximity between individuals.

Formally, we define a network as a graph $G=(V,E)$ , where $V=\{1,2,\dots ,N\}$ is the set of nodes, $E=\left\{\right(u,v):u,v\in V\}$ is the set of edges, and each edge $(u,v)$ defines a pair of nodes $u,v\in V$ that are directly connected. The pairs $(u,v)$ may be ordered, in which case edges are directed (by convention from $u$ to $v$ ), or unordered, in which case edges are undirected (i.e., $v$ connects to $u$ if and only if $u$ connects to $v$ ). We may also add weights to the edges to represent different strengths of interactions, or we may even consider graphs that evolve in time (so-called temporal networks) to account for the dynamics of interactions in a population.

We will primarily consider simple (i.e., static, undirected, and unweighted) networks in our examples. However, the Gillespie algorithms apply to simulated jump processes in all kinds of populations and networks. (For temporal networks, we need to extend the classic Gillespie algorithms to cope with the time-varying network structure; see Section 5.4.)

2.3 Example: Stochastic SIR Model in Continuous Time

We introduce jump processes and explore their mathematical properties by way of a running example. We show how we can use them to model epidemic dynamics using the stochastic susceptible-infectious-recovered (SIR) model.Footnote ¹ For more examples (namely, SIR epidemic dynamics in metapopulation networks, the voter model, and the Lotka–Volterra model for predator–prey dynamics), see Section 3.4.

We examine a stochastic version of the SIR model in continuous time defined as follows. We consider a constant population of $N$ individuals (nodes). At any time, each individual is in one of three states: susceptible (denoted by $S$ ; meaning healthy), infectious (denoted by $I$ ), or recovered (denoted by $R$ ). The rules governing how individuals change their states are shown schematically in Fig. 1. An infectious individual that is in contact with a susceptible individual infects the susceptible individual in a stochastic manner with a constant infection rate $\beta$ . Independently of the infection events, an infectious individual may recover at any point in time, with a constant recovery rate $\mu$ . If an infection event occurs, the susceptible individual that has been infected changes its state to I. If an infectious individual recovers, it transits from the I state to the R state. Nobody leaves or joins the population over the course of the dynamics. After reaching the R state, an individual cannot be reinfected or infect others again. Therefore, R individuals do not influence the increase or decrease in the number of S or I individuals. Because R individuals are as if they no longer exist in the system, the R state is mathematically equivalent to having died of the infection; once dead, an individual will not be reinfected or infect others.

Figure 1 Rules of state changes in the SIR model. An infectious individual infects a susceptible neighbor at a rate $\beta$ . Each infectious individual recovers at a rate $\mu$ .

We typically start the stochastic SIR dynamics with a single infectious individual, which we refer to as the source or seed, and $N_{\text{S}}=N-1$ susceptible individuals (and thus no recovered individuals). Then, various infection and recovery events may occur. The dynamics stop when no infectious individuals are left. In this final situation, the population is composed entirely of susceptible and/or recovered individuals. Since both infection and recovery involve an infectious individual, and there are no infectious individuals left, the dynamics are stuck. The final number of recovered nodes, denoted by $N_{\text{R}}$ , is called the epidemic size, also known as the final epidemic size or simply the final size.Footnote ² The epidemic size tends to increase as the infection rate $\beta$ increases or as the recovery rate $\mu$ decreases. Many other measures to quantify the behavior of the SIR model exist (Reference Pastor-Satorras, Castellano, Van Mieghem and VespignaniPastor-Satorras et al., 2015). For example, we may be interested in the time until the dynamics terminate or in the speed at which the number of infectious individuals grows in the initial stage of the dynamics.

Consider Fig. 2(a), where individuals are connected as a network. We generally assume that infection may only occur between pairs of individuals that are directly connected by an edge (called adjacent nodes). For example, the node $v_4$ can infect $v_1$ and $v_5$ but not $v_3$ . The network version of the SIR model is fully described by the infection rate $\beta$ , the recovery rate $\mu$ , the network structure, that is, which node pairs are connected by an edge, and the choice of source node to initialize the dynamics.

Figure 2 Stochastic SIR process on a square-grid network with six nodes. (a) Status of the network at an arbitrary time $t$ . (b) Status of the network after $v_4$ has recovered. The values attached to the nodes indicate the rates of the events that the nodes may experience next.

Mathematically, we describe the system by a set of coupled, constant-rate jump processes; constant-rate jump processes are known as Poisson processes (Box 1). Each possible event that may happen is associated to a Poisson process, that is, the recovery of each infectious individual is described by a Poisson process, and so is each pair of infectious and susceptible individuals where the former may infect the latter. The Poisson processes are coupled because an event generated by one process may alter the other processes by changing their rates, generating new Poisson processes, or making existing ones disappear. For example, after a node gets infected, it may in turn infect any of its susceptible neighbors, which we represent mathematically by adding new Poisson processes. This coupling implies that the set of coupled Poisson processes generally constitutes a process that is more complicated than a single Poisson process.

In the following subsections we develop the main mathematical properties of Poisson processes and of sets of Poisson processes. We will rely on these properties in Section 3 to construct the Gillespie algorithms that can simulate systems of coupled Poisson processes exactly. Note that the restriction to Poisson (i.e., constant-rate) processes is essential for the classic Gillespie algorithms to work; see Section 5 for recent extensions to the simulation of non-Poissonian processes.

Box 1 Properties of Poisson Processes

A Poisson process is a jump process that generates events with a constant rate, $\lambda$ .

Waiting-Time Distribution

The waiting times $\tau$ between consecutive events generated by a Poisson process are exponentially distributed. In other words, $\tau$ obeys the probability density

$\psi \left(\tau \right)=\lambda e^{-\lambda \tau }.$(2.1)

Memoryless Property

The waiting time left until a Poisson process generates an event given that a time $t$ has already elapsed since the last event is independent of $t$ . This property is called the memoryless property of Poisson processes and is shown as follows:

$\psi (t+\tau |t)=\frac{\psi (t+\tau )}{\Psi \left(t\right)}=\frac{\lambda e^{-\lambda (t+\tau )}}{e^{-\lambda t}}=\lambda e^{-\lambda \tau },$(2.2)

where $\psi (t+\tau |t)$ represents the conditional probability density that the next event occurs a time $t+\tau$ after the last event given that time $t$ has already elapsed; $\Psi \left(t\right)={\int }_t^{\infty }\psi \left(\tau \right)\text{d}\tau =e^{-\lambda t}$ is called the survival probability and is the probability that no event takes place for a time $t$ . The first equality in Eq. (2.2) follows from the definition of the conditional probability. The second equality follows from Eq. (2.1).

Superposition Theorem

Consider a set of Poisson processes indexed by $i\in \{1,2,\dots ,M\}$ . The superposition of the processes is a jump process that generates an event whenever any of the individual processes does. It is another Poisson process whose rate is given by

$\Lambda =\sum _{i=1}^M{\lambda }_i,$(2.3)

where ${\lambda }_i$ is the rate of the $i$ th Poisson process.

Probability of a Given Process Generating an Event in a Superposition of Poisson Processes

Consider any given event generated by a superposition of Poisson processes. The probability ${\Pi }_i$ that the $i$ th individual Poisson process has generated this event is proportional to the rate of the $i$ th process. In other words,

${\Pi }_i={\lambda }_i/\Lambda .$(2.4)

2.4 Waiting-Time Distribution for a Poisson Process

We derive in this subsection the waiting-time distribution for a Poisson process, which characterizes how long one has to wait for the process to generate an event. It is often easiest to start from a discrete-time description when exploring properties of a continuous-time stochastic process. Therefore, we will follow this approach here. We use the recovery of a single node in the SIR model as an example in our development.

Let us partition time into short intervals of length $\delta t$ . As $\delta t$ goes to zero, this becomes an exact description of the continuous-time process. An infectious individual recovers with probability $\mu \delta t$ after each interval given that it has not recovered before.Footnote ³

Formally, we define the SIR process in the limit $\delta t\to 0$ . Then, you might worry that the recovery event is unlikely to ever take place because the probability with which it happens during each time-step, that is, $\mu \delta t$ , goes toward 0 when the step size $\delta t$ does so. However, this is not the case; because the number of time-steps in any given finite interval grows inversely proportional to $\delta t$ , the probability to recover in finite time stays finite. For example, if we use a different step size $\overline{\delta t}=\delta t/10$ , which is ten times smaller than the original $\delta t$ , then the probability of recovery within the short duration of time $\overline{\delta t}$ is indeed 10 times smaller than $\mu \delta t$ (i.e., $=\mu \overline{\delta t}$ ). However, there are $\delta t/\overline{\delta t}=10$ windows of size $\overline{\delta t}$ in one time window of size $\delta t$ . So, we now have 10 chances for recovery to happen instead of only one chance. The probability for recovery to occur in any of these 10 time windows is equal to one minus the probability that it does not occur. The probability that the individual does not recover in time $\delta t$ is equal to $(1-\mu \overline{\delta t}{)}^{\delta t/\overline{\delta t}}$ . Therefore, the probability that the individual recovers in any of the $\delta t/\overline{\delta t}$ windows is

$p_{I\to R}=1-(1-\mu \overline{\delta t}{)}^{\delta t/\overline{\delta t}}.$(2.5)

Equation (2.5) does not vanish as we make $\overline{\delta t}$ small. In fact, the Taylor expansion of Eq. (2.5) in terms of $\overline{\delta t}$ yields $p_{I\to R}\approx \left(\delta t/\overline{\delta t}\right)\times \mu \overline{\delta t}=\mu \delta t$ , where $\approx$ represents “approximately equal to”. Therefore, to leading order, the recovery probabilities are the same between the case of a single time window of size $\delta t$ and the case of $\delta t/\overline{\delta t}$ time windows of size $\overline{\delta t}$ .

In the limit $\delta t\to 0$ , the recovery event may happen at any continuous point in time. We denote by $\tau$ the waiting time from the present time until the time of the recovery event. We want to determine the probability density function (probability density or pdf for short) of $\tau$ , which we denote by ${\psi }_{\text{I}\to \text{R}}\left(\tau \right)$ . By definition, ${\psi }_{\text{I}\to \text{R}}\left(\tau \right)\delta t$ is equal to the probability that the recovery event happens in the interval $[\tau ,\tau +\delta t)$ for an infinitesimal $\delta t$ (i.e., for $\delta t\to 0$ ). To calculate ${\psi }_{\text{I}\to \text{R}}\left(\tau \right)$ , we note that the probability that the event occurs after $r=\tau /\delta t$ time windows, denoted by $p_{\text{I}\to \text{R}}\left(r\right)$ , is equal to the probability that it did not occur during the first $r$ time windows and then occurs in the $(r+1)$ th window. This probability is equal to

$p_{\text{I}\to \text{R}}\left(r\right)=(1-\mu \delta t{)}^r\times \mu \delta t=(1-\mu \delta t{)}^{\tau /\delta t}\mu \delta t.$(2.6)

The first factor on the right-hand side of Eq. (2.6) is the probability that the event has not happened before the $(r+1)$ th window; it is simply equal to the probability that the event has not happened during a single window, raised to the power of $r$ . The second factor is the probability that the event happens in the $(r+1)$ th window. By applying the identity , known from calculus (see Appendix), with $x=-\mu \delta t$ to Eq. (2.6), we obtain the pdf of the waiting time as follows:

$\begin{array}{cc}{\psi }_{\text{I}\to \text{R}}\left(\tau \right)& =\underset{\delta t\to 0}{lim}\frac{p_{\text{I}\to \text{R}}\left(\tau /\delta t\right)}{\delta t}\\ & =\mu \underset{\delta t\to 0}{lim}(1-\mu \delta t{)}^{\tau /\delta t}\\ & =\mu {\left(\underset{\delta t\to 0}{lim}(1-\mu \delta t{)}^{1/(-\mu \delta t)}\right)}^{-\mu \tau }\\ & =\mu e^{-\mu \tau }.\end{array}$(2.7)

Equation (2.7) shows the intricate connection between the Poisson process and the exponential distribution: the waiting time of a Poisson process with rate $\mu$ (here, specifically the recovery rate) follows an exponential distribution with rate $\mu$ (Box 1). This fact implies that the mean time we have to wait for the recovery event to happen is $1/\mu$ . The exponential waiting-time distribution actually completely characterizes the Poisson process. In other words, the Poisson process is the only jump process that generates events separated by waiting times that follow a fixed exponential distribution.

If we consider the infection process between a pair of S and I nodes in complete isolation from the other infection and recovery processes in the population, then exactly the same argument (Eq. (2.7)) holds true. In other words, the time until infection takes place between the two nodes is exponentially distributed with rate $\beta$ , that is,

${\psi }_{\text{S}\to \text{I}}\left(\tau \right)=\beta e^{-\beta \tau }.$(2.8)

However, in practice the infection process is more complicated than the recovery process because it is coupled to other processes. Specifically, if another process generates an event before the infection process does, then Eq. (2.8) may no longer hold true for the infection process in question. For example, consider a node $v_1$ that is currently susceptible and an adjacent node $v_2$ that is infectious, as in Fig. 2. For this pair of nodes, two events are possible: $v_2$ may infect $v_1$ , or $v_2$ may recover. As long as neither of the events has yet taken place, either of the two corresponding Poisson processes may generate an event at any point in time, following Eqs. (2.8) and (2.7), respectively. However, if $v_2$ recovers before it infects $v_1$ , then the infection event is no longer possible, and so Eq. (2.8) no longer holds. We explore in the following two subsections how to mathematically deal with this coupling.

2.5 Independence and Interdependence of Jump Processes

Most models based on jump processes and most simulation methods, including the Gillespie algorithms, implicitly assume that different concurrent jump processes are independent of each other in the sense that the internal state of one process does not influence another. This notion of independence may be a source of confusion because a given process may depend on the events generated earlier by other processes, that is, the processes may be coupled, as we saw is the case for the infection processes in the SIR model. In this section, we sort out the notions of independence and coupling and what they mean for the types of jump processes we want to simulate. We will also explore another type of independence of Poisson processes, which is their independence of the past, called the memoryless property.

We can state the independence assumption as the condition that different processes are only allowed to influence each other by changing the state of the system. In other words, at any point in time each process generates an event at a rate that is independent of all other processes given the current state of the system, that is, the processes are conditionally independent. For example, the rate at which $v_2$ infects $v_1$ in Fig. 2(a) depends on $v_2$ being infectious and $v_1$ being susceptible (corresponding to the system’s current state). However, it does not depend on any internal state of $v_2$ ’s recovery process such as the time left until $v_2$ recovers. Given the states of all nodes, the two processes are independent. Poisson processes are always conditionally independent in this sense. The conditional independence property follows directly from the fact that Poisson processes have constant rates by definition and thus are not influenced by other processes. The conditional independence is essential for the Gillespie algorithms to work. Even the recent extensions of the Gillespie algorithms to simulate non-Poissonian processes, which we review in Section 5, rely on an assumption of conditional independence between the jump processes.

We underline that the assumption of conditional independence does not imply that the different jump processes are not coupled with each other. Such uncoupled processes would indeed be boring. If the jump processes constituting a given system were all uncoupled, then they would not be able to generate any collective dynamics. On the technical side, there would in this case be no reason to consider the set of processes as one system. It would suffice to analyze each process separately. We would in particular have no need for the specialized machinery of the Gillespie algorithms since we could simply simulate each process by sampling waiting times from the corresponding exponential distribution (Box 1, Eq. (2.1)).

In fact, the conditional independence assumption allows different processes to be coupled, as long as they only do so by changing the physical state of the system. This is a natural constraint in many systems. For example, in chemical reaction systems, the processes (that is, chemical reactions) are coupled through discrete reaction events that use molecules of some chemical species to generate others. Similarly, in the SIR model different processes influence each other by changing the state of the nodes, that is, from S to I in an infection event or from I to R in a recovery event. In the example shown in Fig. 2, when node $v_4$ recovers, it decreases the probability that its neighboring susceptible node $v_1$ gets infected within a certain time horizon compared with the scenario where $v_4$ remains infectious. As this example suggests, the probability that a susceptible node gets infected depends on the past states of its neighbors. Therefore, over the course of the entire simulation, the dynamics of a node’s state (e.g., $v_1$ ) are dependent on those of its neighbors (e.g., $v_2$ and $v_4$ ).

Because of the coupling between jump processes, which is present in most systems of interest, we cannot simply simulate the system by separately generating the waiting times for each process according to Eq. (2.1). Any event that occurs will alter the processes to which it is coupled, thus rendering the waiting times we drew for the affected processes invalid. What the Gillespie algorithms do instead is to successively generate the waiting time until the next event, update the state of the system, and reiterate.

Besides being conditionally independent of each other, Poisson processes also display a temporal independence property, the so-called memoryless property (Box 1). In Poisson processes, the probability of the time to the next event, $\tau$ , is independent of how long we have already waited since the last event. In this sense, we do not need to worry about what has happened in the past. The only things that matter are the present status of the population (such as $v_1$ is susceptible and $v_2$ is infectious right now) and the model parameters (such as $\beta$ and $\mu$ ). The memoryless property can be seen as a direct consequence of the exponential distribution of waiting times of Poisson processes (Box 1, Eq.(2.2)).

2.6 Superposition of Poisson Processes

In this section, we explain a remarkable property of Poisson processes called the superposition theorem. The direct method exploits this theorem. Other methods, such as the rejection sampling algorithm (see Section 2.8) and the first reaction method, can also benefit from the superposition theorem to accelerate the simulations without impacting their accuracy.

Consider a susceptible individual $v_i$ in the SIR model that is in contact with $N_{\text{I}}$ infectious individuals. Any of the $N_{\text{I}}$ infectious individuals may infect $v_i$ . Consider the case shown in Fig. 3(a), where $N_{\text{I}}=3$ . If we focus on a single edge connecting $v_i$ to one of its neighbors and ignore the other neighbors, the probability that $v_i$ is infected via this edge exactly in time $[\tau ,\tau +\delta t)$ from now, where $\delta t$ is small, is given by ${\psi }_{\text{S}\to \text{I}}\left(\tau \right)\delta t$ (see Eq. (2.8)). Each of $v_i$ ’s $N_{\text{I}}$ neighbors may infect $v_i$ in the same manner and independently. The neighbor that does infect $v_i$ is the one for which the corresponding waiting time is the shortest, provided that it does not recover before it infects $v_i$ . From this we can intuitively see that the larger $N_{\text{I}}$ is, the shorter the waiting time before $v_i$ gets infected tends to be. To simulate the dynamics of this small system, we need to know, not when each of its neighbors would infect $v_i$ , but rather the time until any of its neighbors infects $v_i$ .

Figure 3 A susceptible node and other nodes surrounding it. (a) A susceptible node $v_i$ surrounded by three infectious nodes. (b) A susceptible node $v_i$ surrounded by five nodes in different states.

To calculate the waiting-time distribution for the infection of $v_i$ by any of its neighbors, we again resort to the discrete-time view of the infection processes. Because the infection processes are independent, the probability that $v_i$ is not infected by any of its $N_{\text{I}}$ infectious neighbors in a time window of duration $\delta t$ is given by

$(1-\beta \delta t{)}^{N_{\text{I}}}.$(2.9)

Therefore, the probability that $v_i$ is infected after a time $\tau =r\delta t$ (i.e., $v_i$ gets infected exactly in the $(r+1)$ th time window of length $\delta t$ and not before) is given by

$p_{\text{I}\to \text{R}}={\left((1-\beta \delta t{)}^{N_{\text{I}}}\right)}^r\times \left(1-(1-\beta \delta t{)}^{N_{\text{I}}}\right).$(2.10)

Here, the factor ${\left((1-\beta \delta t{)}^{N_{\text{I}}}\right)}^r$ is the survival probability that an infection does not happen for a time $\tau =r\delta t$ . The factor $\left(1-(1-\beta \delta t{)}^{N_{\text{I}}}\right)$ is the probability that any of $v_i$ ’s infectious neighbors infects $v_i$ in the next time window, $t\in [\tau ,\tau +\delta t)$ .

Using the exponential identity $\underset{x\to 0}{lim}(1+x{)}^{1/x}=e$ with $x=-\beta \delta t$ as we did in Section 2.4, we obtain in the continuous-time limit that

$\underset{\delta t\to 0}{lim}{\left((1-\beta \delta t{)}^{N_{\text{I}}}\right)}^r=\underset{\delta t\to 0}{lim}{\left((1-\beta \delta t{)}^{1/(-\beta \delta t)}\right)}^{-N_{\text{I}}\beta \tau }=e^{-N_1\beta \tau },$(2.11)

where the first equality is obtained by noting that $r=\tau /\delta t$ and rearranging the terms. In the same limit of $\delta t\to 0$ , we obtain from Taylor expansion that

$\left(1-(1-\beta \delta t{)}^{N_{\text{I}}}\right)\approx 1-(1-N_{\text{I}}\beta \delta t)=N_{\text{I}}\beta \delta t.$(2.12)

By combining Eqs. (2.10), (2.11), and (2.12), we obtain $p_{\text{I}\to \text{R}}\approx N_{\text{I}}\beta e^{-N_{\text{I}}\beta \tau }\delta t$ . Therefore, the probability density with which $v_i$ gets infected at time $\tau$ is given by

${\psi }_{\text{I}\to \text{R}}\left(\tau \right)=N_{\text{I}}\beta e^{-N_{\text{I}}\beta \tau },$(2.13)

that is, the exponential distribution with rate parameter $N_{\text{I}}\beta$ . By comparing Eqs. (2.8) and (2.13), we see that the effect of having $N_{\text{I}}$ infectious neighbors (see Fig. 3(a) for the case of $N_{\text{I}}=3$ ) is the same as that of having just one infectious neighbor with an infection rate of $N_{\text{I}}\beta$ .

This is a convenient property of Poisson processes, known as the superposition theorem (see Box 1, Eq. (2.3) for the general theorem). To calculate how likely it is that a susceptible node $v_i$ will be infected in time $\tau$ , one does not need to examine when the infection would happen or whether the infection happens for each of the infectious individuals contacting $v_i$ . We are allowed to agglomerate all those effects into one infectious supernode as if the supernode infects $v_i$ with rate $N_{\text{I}}\beta$ . We refer to such a superposed Poisson process that induces a particular state transition in the system (in the present case, the transition from the S to the I state for $v_i$ ) as a reaction channel, following the nomenclature in chemical reaction systems.

This interpretation remains valid even if $v_i$ is adjacent to other irrelevant individuals. In the network shown in Fig. 3(b), the susceptible node $v_i$ has degree (i.e., number of other nodes that are connected to $i$ by an edge) $k_i=5$ . Three neighbors of $v_i$ are infectious, one is susceptible, and one is recovered. In this case, $v_i$ will be infected at a rate of $3\beta$ , the same as in the case of $v_i$ in the network shown in Fig. 3(a).

In both cases, we are replacing three instances of the probability density of the time to the next infection event, each given by $\beta e^{-\beta \tau }$ , by a single probability density $3\beta e^{-3\beta \tau }$ . Representing the three infectious nodes by one infectious supernode, that is, one reaction channel, with three times the infection rate is equivalent to superposing the three Poisson processes into one. Figure 4 illustrates this superposition, showing the putative event times generated by each Poisson process as well as those generated by their superposition. The superposition theorem dictates that the superposition is a Poisson process with a rate of $3\beta$ . This in particular means that we can draw the waiting time $\tau$ until the first of the events generated by all the three Poisson processes happens (shown by the double-headed arrow in Fig. 4) directly from the exponential distribution $\psi \left(\tau \right)=3\beta e^{-3\beta \tau }$ . Note that Poisson processes are defined as generating events indefinitely, and for illustrative purposes we show multiple events in Figure 4. However, in the SIR model only the first event in the superposed process will take place in practice. For example, once the event changes the state of $v_i$ from S to I, the node cannot be infected anymore, and therefore none of the three infection processes can generate any more events.

Figure 4 Superposition of three Poisson processes. The event sequence in the bottom is the superposition of the three event sequences corresponding to each of the three edges connecting $v_i$ to its neighbors $v_j$ , $v_{j^{\prime }}$ , and $v_{j^{\prime \prime }}$ . The superposed event sequence generates an event whenever one of the three individual processes does. Note that edges ( $v_i$ , $v_j$ ), ( $v_i$ , $v_{j^{\prime }}$ ), and ( $v_i$ , $v_{j^{\prime \prime }}$ ) generally carry different numbers of events in a given time window despite the rate of the processes (i.e., the infection rate, $\beta$ ) being the same. This is due to the stochastic nature of Poisson processes.

Let us consider again the snapshot of the SIR dynamics shown in Fig. 3(a), but this time we consider all the possible infection and recovery events. We can represent all the possible events that may occur by four reaction channels (i.e., Poisson processes). One channel represents the infection of the node $v_i$ by any of its neighbors, which happens at a rate $3\beta$ . We refer to this reaction channel as the first reaction channel. The three other channels each represent the recovery process of one of the infectious nodes. We refer to these three reaction channels as the second to the fourth reaction channels. We can use the same approach as above to obtain the probability density for the waiting time until the first event generated by any of the channels. However, to completely describe the dynamics, it is not sufficient to know when the next event happens. We also need to know which channel generates the event. Precisely speaking, we need to know the probability ${\Pi }_i$ that it is the $i$ th reaction channel that generates the event. Using the definition of conditional probability, we obtain

$\begin{array}{cc}& {\Pi }_i=\frac{\left\{\text{probability that an event in the }i\text{th reaction channel occurs}\right\}}{\left\{\text{probability that an event in any reaction channel }j\in \right\{1,2,3,4\left\}\text{ occurs}\right\}}.\end{array}$(2.14)

In a discrete-time description, the numerator in Eq. (2.14) is simply ${\lambda }_i\delta t$ , where ${\lambda }_1=3\beta$ and ${\lambda }_2={\lambda }_3={\lambda }_4=\mu$ are the rates of the reaction channels. The denominator is equal to $1-\prod _{j=1}^4(1-{\lambda }_j\delta t)$ , which in the limit of small $\delta t$ can be Taylor expanded to $\sum _{j=1}^4{\lambda }_j\delta t=3(\beta +\mu )\delta t$ . Thus, the probability that the $i$ th reaction channel has generated an event that has taken place is

${\Pi }_i=\frac{{\lambda }_i}{\sum _{j=1}^4{\lambda }_j},$(2.15)

that is, ${\Pi }_i$ is simply proportional to the rate ${\lambda }_i$ .

The same result holds true for general superpositions of Poisson processes (see Box 1).

2.7 Ignoring Stochasticity: Differential Equation Approach

We have introduced the types of models we are interested in and have explored their basic mathematical properties. We now turn our attention to the problem of how we can solve such models in practice. We consider again the SIR model. One simple strategy to solve it is to forget about the true stochastic nature of infection and recovery and approximate the processes as being deterministic. In this approach, we only track the dynamics of the mean numbers of susceptible, infectious, and recovered individuals. Such deterministic dynamics are described by a system of ordinary differential equations (ODEs). The ODE version of the SIR model has a longer history than the stochastic one, dating back to the seminal work by William Ogilvy Kermack and Anderson Gray McKendrick in the 1920s (Reference Kermack and McKendrickKermack & McKendrick, 1927). For the basic SIR model described previously, the corresponding ODEs are given by

$\begin{array}{cc}\frac{\text{d}{\rho }_{\text{S}}}{\text{d}t}=& -\beta {\rho }_{\text{S}}{\rho }_{\text{I}},\end{array}$(2.16)

$\begin{array}{cc}\frac{\text{d}{\rho }_{\text{I}}}{\text{d}t}=& \beta {\rho }_{\text{S}}{\rho }_{\text{I}}-\mu {\rho }_{\text{I}},\end{array}$(2.17)

$\begin{array}{cc}\frac{\text{d}{\rho }_{\text{R}}}{\text{d}t}=& \mu {\rho }_{\text{I}},\end{array}$(2.18)

where ${\rho }_{\text{S}}=N_S/N$ , ${\rho }_{\text{I}}=N_I/N$ , and ${\rho }_{\text{R}}=N_R/N$ are the fractions of S, I, and R individuals, respectively. The $\beta {\rho }_{\text{S}}{\rho }_{\text{I}}$ terms in Eqs. (2.16) and (2.17) represent infection events, through which the number of S individuals decreases and the number of I individuals increases by the same amount. The $\mu {\rho }_{\text{I}}$ terms in Eqs. (2.17) and (2.18) represent recovery events.

One can solve Eqs. (2.16), (2.17), and (2.18) either analytically, to some extent, or numerically using an ODE solver implemented in various programming languages. Suppose that we have coded up Eqs. (2.16), (2.17), and (2.18) into an ODE solver to simulate the infection dynamics (such as time courses of ${\rho }_{\text{I}}$ ) for various values of $\beta$ and $\mu$ . Does the result give us complete understanding of the original stochastic SIR model? The answer is negative (Reference Mollison, Isham and GrenfellMollison, Isham, & Grenfell, 1994), at least for the following reasons.

First, the ODE is not a good approximation when $N$ is small. In Eqs. (2.16), (2.17), and (2.18), the variables are the fraction of individuals in each state. For example, ${\rho }_{\text{I}}=N_{\text{I}}/N$ . The ODE description assumes that ${\rho }_{\text{I}}$ can take any real value between 0 and 1 and that ${\rho }_{\text{I}}$ changes continuously as time goes by. However, in reality ${\rho }_{\text{I}}$ is quantized, so it can take only the values $0$ , $1/N$ , $2/N$ , $\dots$ $(N-1)/N$ , and $1$ , and it changes in steps of $1/N$ (e.g., it changes from $3/N$ to $4/N$ discontinuously). This discrete nature does not typically cause serious problems when $N$ is large, in which case ${\rho }_{\text{ I}}$ is close to being continuous. By contrast, the ODE model is not accurate when $N$ is small due to the quantization effect. (Note that the ODE approach is problematic in some cases even when $N$ is large, i.e. near critical points, as we discuss in what follows.)

Second, even if $N$ is large, the actual dynamic changes in ${\rho }_{\text{I}}$ , for example, are not close to what the ODEs describe when the number of infectious individuals is small. For example, if ${\rho }_{\text{I}}=2/N$ , there are two infectious individuals. If one of them recovers, ${\rho }_{\text{I}}$ changes to $1/N$ , and this is a 50 percent decrease in ${\rho }_{\text{I}}$ . The ODE assumes that ${\rho }_{\text{ I}}$ changes continuously and is not ready to describe such a change. As another example, suppose that we initially set ${\rho }_{\text{S}}=(N-1)/N$ , ${\rho }_{\text{I}}=1/N$ , and ${\rho }_{\text{R}}=0$ . In other words, there is a single infectious seed, and all the other individuals are initially susceptible. In fact, the theory of the ODE version of the SIR model shows that ${\rho }_{\text{I}}$ increases deterministically, at least initially, if $\beta >\mu$ , corresponding to the situation in which an outbreak of infection happens. However, in the stochastic SIR model, the only initially infectious individual may recover before it infects anybody even if $\beta >\mu$ . When this situation occurs, the dynamics terminate once the initially infectious individual has recovered, and no outbreak is observed. Although the probability with which this situation occurs decreases as $\beta /\mu$ increases, it is still not negligibly small for many large $\beta /\mu$ values. This is inconsistent with the prediction of the ODE model. It should be noted that another common way to initialize the system is to start with a small fraction of infectious individuals, regardless of $N$ . In this case, if we start the stochastic SIR dynamics in a large, well-mixed population and, for example, with 10 percent initially infectious individuals, the ODE version is sufficiently accurate at describing the stochastic SIR dynamics.

Third, ODEs are not accurate at describing the counterpart stochastic dynamics when the system is close to a so-called critical point. For example, in the SIR model, given the value of the infection rate (i.e., $\mu$ ), there is a value of the infection rate called the epidemic threshold, which we denote by ${\beta }_{\text{c}}$ . For $\beta <{\beta }_{\text{c}}$ , only a small number of the individuals will be infected (i.e., the final epidemic size is of $O\left(1\right)$ ). For $\beta >{\beta }_{\text{c}}$ , the final epidemic size is large (i.e., $O\left(N\right)$ ) with a positive probability. In analogy with statistical physics, ${\beta }_{\text{c}}$ is termed a critical point of the SIR model. Near criticality the fluctuations of ${\rho }_S$ , ${\rho }_I$ , and ${\rho }_R$ are not negligible compared to their mean values, even for large $N$ , and the ODE generally fails.

Fourth, ODEs are not accurate when dynamics are mainly driven by stochasticity rather than by the deterministic terms on the right-hand sides of the ODEs. This situation may happen even far from criticality or in a model that does not show critical dynamics. The voter model (see Section 4.10.3 for details) is such a case. In its simplest version, the voter model describes the tug-of-war between two equally strong opinions in a population of individuals. Because the two opinions are equally strong, the ODE version of the voter model predicts that the fraction of individuals supporting opinion A (and that of individuals supporting opinion B) does not vary over time, that is, one obtains $\text{d}{\rho }_{\text{A}}/\text{d}t=\text{d}{\rho }_{\text{B}}/\text{d}t=0$ , where ${\rho }_{\text{A}}$ and ${\rho }_{\text{B}}$ are the fractions of individuals supporting opinions A and B, respectively. However, in fact, the opinion of the individuals flips here and there in the population due to stochasticity, and it either increases or decreases over time.

To summarize, when stochasticity manifests itself, the approximation of the original stochastic dynamics by an ODE model is not accurate.

2.8 Rejection Sampling Algorithm

The most intuitive method to simulate the stochastic SIR model, while accounting for the stochastic nature of the model, is probably to discretize time and simulate the dynamics by testing whether each possible event takes place in each step. This is called the rejection sampling algorithm. Let us consider the stochastic SIR model on a small network composed of $N=6$ nodes, as shown in Fig. 2, to explain the procedure.

Assume that the state of the network (i.e., the states of the individual nodes) is as shown in Fig. 2(a) at time $t$ ; three nodes are susceptible, two nodes are infectious, and the other node is recovered. In the next time-step, which accounts for a time length of $\Delta t$ and corresponds to the time interval $[t,t+\Delta t)$ , an infection event may happen in five ways: $v_2$ infects $v_1$ , $v_2$ infects $v_3$ , $v_2$ infects $v_5$ , $v_4$ infects $v_1$ , and $v_4$ infects $v_5$ . Recovery events may happen for $v_2$ and $v_4$ . Therefore, there are seven possible events in total, some of which may simultaneously happen in the next time-step.

With the rejection sampling method, we sequentially (called asynchronous updating) or simultaneously (called synchronous updating) check whether or not each of these events happens in each time-step of length $\Delta t$ . Note that it is not possible to go to the limit of $\Delta t\to 0$ in rejection sampling. In our example, $v_2$ infects $v_1$ with probability $\beta \Delta t$ in a time-step. With probability $1-\beta \Delta t$ , nothing occurs along this edge. In practice, to determine whether the event takes place or not, we draw a random number $u$ uniformly from $[0,1)$ . If $u\ge \beta \Delta t$ , the algorithm rejects the proposed infection event (thus the name rejection sampling). If $u<\beta \Delta t$ , we let the infection occur. Then, under asynchronous updating, we change the state of $v_1$ from S to I and update the set of possible events accordingly right away, and then proceed to check the occurrence of each of the remaining possible events in turn. Under synchronous updating, we first check whether each of the possible state changes takes place and note down the changes that take place. We then implement all the noted changes simultaneously. Regardless of whether we use asynchronous or synchronous updating, the infection event occurs with probability $\beta \Delta t$ .

If $v_4$ recovers, which occurs with probability $\mu \Delta t$ , and none of the other six possible events occurs in the same time-step, the status of the network at time $t+\Delta t$ is as shown in Fig. 2(b). Then, in the next time-step, $v_1$ may get infected, $v_2$ may recover, $v_3$ may get infected, and $v_5$ may get infected, which occurs with probabilities $\beta \Delta t$ , $\mu \Delta t$ , $\beta \Delta t$ , and $\beta \Delta t$ , respectively. In this manner, we carry forward the simulation by discrete steps until no infectious nodes are left.

There are several caveats to this approach. First, the asynchronous and the synchronous updating schemes of the same stochastic dynamics model may lead to systematically different results (Reference Cornforth, Green and NewthCornforth, Green, & Newth, 2005; Reference Greil and DrosselGreil & Drossel, 2005; Reference Huberman and GlanceHuberman & Glance, 1993).

Second, one should set $\Delta t$ such that both $\beta \Delta t<1$ and $\mu \Delta t<1$ always hold true. In fact, the discrete-time interpretation of the original model is justified only when $\Delta t$ is small enough to yield $\beta \Delta t\ll 1$ and $\mu \Delta t\ll 1$ .

Third, in the case of asynchronous updating, the order of checking the events is arbitrary, but it affects the outcome, particularly if $\Delta t$ is not tiny. For example, we can sequentially check whether each of the five infection events occurs and then whether each of the two recovery events occurs, completing one time-step, One can alternatively check the recovery events first and then the infection events. If we do so and $v_4$ recovers in the time-step, then it is no longer possible that $v_4$ infects $v_1$ or $v_5$ in the same time-step because $v_4$ has recovered. If the infection events were checked before the recovery events, it is possible that $v_4$ infects $v_1$ or $v_5$ before $v_4$ recovers in the same time-step.

Fourth, some of the seven types of event cannot occur simultaneously in a single time-step regardless of whether the updating is asynchronous or synchronous, and regardless of the order in which we check the events in the asynchronous updating. For example, if $v_2$ has infected $v_1$ , then $v_4$ cannot infect $v_1$ in the same time-step (or anytime later) and vice versa. In fact, from the susceptible node $v_1$ ’s point of view, it does not matter which infectious neighbor, either $v_2$ or $v_4$ , infects $v_1$ . What is primarily important is whether $v_1$ gets infected or not in the given time-step, whereas one wants to know who infected whom in some tasks such as contact tracing.

A useful method to mitigate the effect of overlapping events of this type is to take a node-centric view. The superposition theorem implies that $v_1$ will get infected according to a Poisson process with rate $2\beta$ because it has two infectious neighbors (Section 2.6). By exploiting this observation, let us redefine the list of possible events at time $t$ . The node $v_1$ will get infected with probability $2\beta \Delta t$ (and will not get infected with probability $1-2\beta \Delta t$ ). Nodes $v_3$ and $v_5$ will get infected with probabilities $\beta \Delta t$ and $2\beta \Delta t$ , respectively. As before, $v_2$ and $v_4$ recover with probability $\mu \Delta t$ each. In this manner, we have reduced the number of possible events from seven to five. We are usually interested in simulating such stochastic processes in much larger networks or populations, where nodes tend to have a degree larger than in the network shown in Fig. 2. For example, if a node $v_i$ has 50 infectious neighbors, implementing the rejection sampling using the probability that $v_i$ gets infected, $50\beta \Delta t$ , rather than checking if $v_i$ gets infected with probability $\beta \Delta t$ along each edge that $v_i$ has with an infectious neighbor, will confer a fiftyfold speedup of the algorithm.

Rejection sampling is a widely used method, particularly in research communities where continuous-time stochastic process thinking does not prevail. In a related vein, many people are confused by being told that the infection and recovery rates $\beta$ and $\mu$ can exceed 1. They are accustomed to think in discrete time such that they are not trained to distinguish between the rate and probability. They are different; simply put, the rate is for continuous time, and the probability is for discrete time. Here we advocate that we should not use the discrete-time versions in general, despite their simplicity and their appeal to our intuition, for the following reasons (see Reference Gómez, Gómez-Gardeñes, Moreno and ArenasGómez, Gómez-Gardeñes, Moreno, and Arenas [2011] and Reference Fennell, Melnik and GleesonFennell, Melnik, and Gleeson [2016] for similar arguments).

First, the use of a small $\Delta t$ , which is necessary to assure an accurate approximation of the actual continuous-time stochastic process, implies a large computation time. If the duration of time that one run of simulation needs is $T$ , one needs $n=T/\Delta t$ discrete time-steps, which is large when $\Delta t$ is small. How small should $\Delta t$ be? It is difficult to say. If you run simulations with a choice of a small $\Delta t$ and calculate statistics of your interest or draw a figure for your report or paper, a good practice is to try the same thing after halving $\Delta t$ . If the results do not noticeably change, then your original choice of $\Delta t$ is probably small enough for your purpose. Otherwise, you need to make $\Delta t$ smaller. It takes time to carry out such a check just to determine an appropriate $\Delta t$ value. Many people skip it. The Gillespie algorithms do not rely on a discrete-time approximation and are also typically faster than rejection sampling with a reasonably small $\Delta t$ value.

Second, no matter how small $\Delta t$ is, the results of rejection sampling are only approximate. This is because it is exact only in the limit $\Delta t\to 0$ . By contrast, the Gillespie algorithms are always exact.

Proponents of the rejection sampling method may say that they want to define the model (such as the SIR model) in discrete time and run it rather than to consider the continuous-time version of the model and worry about the choice of $\Delta t$ or the accuracy of the rejection sampling. We recommend against this as well. In the SIR model in discrete time, any infectious individual $v_i$ infects a neighboring susceptible individual $v_j$ with probability ${\beta }^{\prime }$ , and each infectious individual recovers with probability ${\mu }^{\prime }$ . Then, there are at least two problems related to this. First, the order of the events affects dynamics in the case of asynchronous updating. Second, and more importantly, we do not know how to change the time resolution of the simulation when we need to. For example, if one simulation step currently corresponds to one hour, one may want to now simulate the same model with some more temporal detail such that one step corresponds to ten minutes. Because the physical time is now one sixth of the original one, should we multiply ${\beta }^{\prime }$ and ${\mu }^{\prime }$ by $1/6$ and do the same simulations? The answer is no. If the original time-step corresponds to $\Delta t=1$ , which is often implicit, then the probability that an infectious individual recovers in the continuous-time stochastic SIR model within time $\Delta t(=1)$ is $1-e^{-\mu \Delta t}=1-e^{-\mu }$ , which we equate with ${\mu }^{\prime }$ . Then, if we scale the time $c$ times (e.g., $c=1/6$ ), the probability that the recovery occurs in a new single time-step is $1-e^{-\mu c\Delta t}=1-e^{-c\mu }$ , which is not equal to $c{\mu }^{\prime }$ . For example, with $\mu =1$ and $c=1/6$ , one obtains $1-e^{-c\mu }\approx 0.154$ , whereas $c{\mu }^{\prime }\approx 0.105$ .

3.1 Brief History

As the name suggests, the Gillespie algorithms are ascribed to American physicist Daniel Thomas Gillespie (Reference GillespieGillespie, 1976, Reference Gillespie1977). He originally proposed them in 1976 for simulating stochastic chemical reaction systems, and they have seen many applications as well as further algorithmic developments in this field. Nevertheless, the algorithms only rely on general properties of Poisson processes and not on any particular properties of chemical reactions. Therefore, the applicability of the Gillespie algorithms is much wider than to chemical reaction systems. In fact, they have been extensively used in simulations of multiagent systems both in unstructured populations and on networks. The only assumptions are that the system undergoes changes via sequences of discrete events (e.g., somebody infects somebody, somebody changes its internal state from a low-activity state to a high-activity state) and with a rate that stays constant in between events. Nevertheless, the latter assumption has been relaxed in recent extensions of the algorithms; we review these in Section 5.

There were precursors to the Gillespie algorithms. The American mathematician Joseph Leo Doob developed in his 1942 and 1945 papers the mathematical foundations of continuous-time Markov chains that underlie the Gillespie algorithms (Reference DoobDoob, 1942, Reference Doob1945). In the second of the two papers, he effectively proposed the direct method, although the focus of the paper was mathematical theory and he did not propose a computational implementation Reference DoobDoob (1945: pp. 465–466). Due to this, the algorithm is sometimes called the Doob–Gillespie algorithm. David George Kendall, who is famous for Kendall’s notation in queuing theory,Footnote ⁴ implemented an equivalent of the direct method to simulate a stochastic birth-death process on a computer as early as in 1950 (Reference KendallKendall, 1950). In 1953, Maurice Stevenson Bartlett, a British statistician, simulated the SIR model in a well-mixed population (i.e., every pair of individuals is directly connected to each other) using the direct method (Reference BartlettBartlett, 1953).

Independently of Gillespie, Alfred B. Bortz and colleagues also proposed the same algorithm as the direct method to simulate stochastic dynamics of Ising spin systems in statistical physics in 1975 (Reference Bortz, Kalos and LebowitzBortz, Kalos, & Lebowitz, 1975). Therefore, the direct method is also called the Bortz–Kalos–Lebowitz algorithm (or the n-fold way following the naming in their paper, and also rejection-free kinetic Monte Carlo and the residence-time algorithm). An even earlier paper published in 1966 in the same field proposed almost the same algorithm, with the only difference being that the waiting time between events was assumed to take a deterministic value rather than being stochastic as in the Gillespie algorithms (Reference Young and ElcockYoung & Elcock, 1966).

3.2 First reaction method

To introduce the first reaction method, we consider our earlier example of the SIR model on a six-node network (see Fig. 2(a)). Here two types of events may happen next: either a susceptible node becomes infected (S $\to$ I), or an infectious node recovers (I $\to$ R). The rates at which each node experiences a state transition are shown in Fig. 5(a), which replicates Fig. 2(a). For example, $v_1$ is twice as likely to be infected next as $v_3$ is because $v_1$ has two infectious neighbors whereas $v_3$ has one infectious neighbor. Each event obeys a separate Poisson process. Therefore, let us first generate hypothetical event sequences according to each Poisson process with their respective rates (see Fig. 5(b)). In fact, we need to use at most only the first event in each sequence (shown in magenta in Fig. 5(b)). For example, in the event sequence for $v_1$ , the first event may be used, in which case $v_1$ will be infected. Once $v_1$ is infected, the subsequent events on the same event sequence will be discarded because $v_1$ will never be infected again. If $v_1$ is infected, it will undergo another type of event, which is recovery. However, we cannot reuse the second or any subsequent events in the same sequence as the recovery event because the recovery occurs with rate $\mu$ , which is different from the rate $2\beta$ with which we have generated the event sequence for $v_1$ . A lesson we learn from this example is that we should not prepare many possible event times beforehand because most of them would be discarded.

Figure 5 Determination of the time to the next event in the SIR model using the Gillespie algorithms. (a) The current state of the system. The rate with which each node changes its state is shown next to the node. (b) The putative events generated by the Poisson process corresponding to each node, with the first event of each process shown in magenta. There is no event on the timeline of $v_6$ . This is because $v_6$ is in the recovered state and no longer undergoes any state change; therefore the event rate of the Poisson process associated with $v_6$ is equal to 0. The first reaction method generates a putative waiting time until the first event for each process and selects the smallest one. The direct method directly generates the waiting time until the first event in the superposed process (bottom of panel (b)).

Given this reasoning, Gillespie proposed the first reaction method in one of his two original papers on the Gillespie algorithms (Reference GillespieGillespie, 1976). The idea of the first reaction method is to generate only the first putative event time (shown in magenta in Fig. 5(b)) for each node. These putative times correspond to when each node would experience a next event if nothing else happens in the rest of the system. We do not generate an event time for $v_6$ because this node is recovered, so it will never undergo any event. We then figure out which event occurs first. (In the example in Fig. 5(a), it is node $v_4$ that will change its state, and the state change is from I to R.)

To generate a putative waiting time for each node $v_1$ , $\dots$ , $v_5$ in practice, we use a general technique called inverse sampling, which proceeds as follows. For example, the time to the first event for $v_1$ obeys the exponential distribution ${\psi }_1\left(\tau \right)\equiv 2\beta e^{-2\beta \tau }$ . The probability that the time to the next event is larger than $\tau$ , called the survival probability (also called the survival function and the complementary cumulative distribution function), is given by

${\Psi }_1\left(\tau \right)\equiv {\int }_{\tau }^{\infty }{\psi }_1\left({\tau }^{\prime }\right)\text{d}{\tau }^{\prime }={\int }_{\tau }^{\infty }2\beta e^{-2\beta {\tau }^{\prime }}\text{d}{\tau }^{\prime }=e^{-2\beta \tau }.$(3.1)

By definition, ${\Psi }_1\left(\tau \right)$ is a probability, so $0<{\Psi }_1\left(\tau \right)\le 1$ . We have excluded ${\Psi }_1\left(\tau \right)=0$ because it happens only in the limit $\tau \to \infty$ . This exclusion does not cause any problem in the following development.

Then, we draw a number $u$ from an unbiased random number generator that generates random numbers uniformly in the interval $(0,1]$ . The generated number $u$ is called a uniform $(0,1]$ random variate. Any practical programming language has a function to generate uniform (pseudo) random variates. We do, however, advise against using a programming language’s standard pseudo-random number generator, which is typically of poor quality. You should instead use one from a scientific programming library or code it yourself. We will discuss some good practices for pseudorandom number generation in Section 4.9. Given $u$ , we then generate $\tau$ from the implicit equation

$u={\Psi }_1\left(\tau \right).$(3.2)

For example, if we draw $u=0.3$ , we can find the unique $\tau$ value that satisfies Eq. (3.2). This method for generating random variates obeying a given distribution is called inverse (transform) sampling (Reference von Neumannvon Neumann, 1951). One can use this method for a general probability density function, $\psi \left(\tau \right)$ , as long as one can calculate its survival function, $\Psi \left(\tau \right)={\int }_{\tau }^{\infty }\psi \left({\tau }^{\prime }\right)\text{d}{\tau }^{\prime }$ , like in Eq. (3.1). In the present case, by combining Eqs. (3.1) and (3.2), we obtain

$u=e^{-2\beta \tau },$(3.3)

which leads to

$\tau =-\frac{lnu}{2\beta }.$(3.4)

Note that $\tau$ $\ge 0$ because $lnu\le 0$ . Equation (3.4) is reasonable in the sense that a large event rate $2\beta$ will yield a small waiting time $\tau$ on average.

In the same manner, one can generate the putative times to the next event for $v_1$ , $\dots$ , $v_5$ , denoted by ${\tau }_1^{\text{put}}$ , $\dots$ , ${\tau }_5^{\text{put}}$ , using five uniform random variates. If the realized ${\tau }_1^{\text{put}}$ , $\dots$ , ${\tau }_5^{\text{put}}$ values are as shown in Fig. 5(b), we conclude that node $v_4$ recovers next. Then, we change the state of $v_4$ from I to R and advance the clock by time $\tau ={\tau }_4^{\text{put}}$ . Once $v_4$ recovers, the configuration of the six nodes will be the one given in Fig. 2(b). We then repeat the same procedure to find the next event given the updated set of processes corresponding to the nodes’ states after the event (i.e., we now have three infection processes with rate $\beta$ and one recovery process with rate $\mu$ ), and so on.

Note that in our example, the event rate changed for $v_1$ , $v_4$ , and $v_5$ , while it remained unchanged for $v_2$ and $v_3$ . Therefore, we do not need to generate entirely new putative waiting times for $v_2$ and $v_3$ . We just have to update ${\tau }_2^{\text{put}}$ and ${\tau }_3^{\text{put}}$ as ${\tau }_2^{\text{put}}\to {\tau }_2^{\text{put}}-\tau$ and ${\tau }_3^{\text{put}}\to {\tau }_3^{\text{put}}-\tau$ , respectively, to account for the time $\tau$ that has elapsed.

In this manner, we can reuse ${\tau }_2^{\text{put}}$ and ${\tau }_3^{\text{put}}$ (by subtracting $\tau$ ) and avoid having to generate new pseudorandom numbers for redrawing ${\tau }_2^{\text{put}}$ and ${\tau }_3^{\text{ put}}$ . The effect of this frugality becomes important for larger systems where an event generally affects only a small fraction of the processes. The so-called next reaction method (Reference Gibson and BruckGibson & Bruck, 2000) exploits this idea to improve the computational efficiency of the first reaction method. Although the classic first reaction method did not make use of this trick, we include it here because it is simple to implement.

Let us go back to our example. For $v_4$ , we no longer need to generate the time to the next event because $v_4$ is now in the R state. For $v_1$ and $v_5$ , we need to discard ${\tau }_1^{\text{put}}$ and ${\tau }_5^{\text{ put}}$ because they were generated under the assumption that the event rate was $2\beta$ . Now, we need to redraw the time to the next event for the two nodes according to the new distribution ${\psi }_1\left(\tau \right)={\psi }_5\left(\tau \right)=\beta e^{-\beta \tau }$ . For example, we reset ${\tau }_1^{\text{put}}=-ln{u}^{\prime }/\beta$ , where $u^{\prime }$ is a new uniform $(0,1]$ random variate. Although the time $\tau$ has passed to transit from the status of the network shown in Fig. 2(a) to that shown in Fig. 2(b), we do not need to take into account the elapsed time (i.e., $\tau$ ) when generating the new ${\tau }_1^{\text{put}}$ and ${\tau }_5^{\text{put}}$ values. This is due to the memoryless property of Poisson processes (see Box 1), that is, what happened in the past, such as how much time has passed to realize the state transition of $v_4$ , is irrelevant.

The first reaction method in its general form is given in Box 2.

3.3 Direct Method

The direct method exploits the superposition theorem to directly generate the waiting times between successive events in the full system of coupled Poisson processes. For expository purposes, we hypothetically generate an event sequence on each node with the respective rate, although only at most the first event in each sequence will be used. We then superpose the nodal event sequences into one Poisson process (see Fig. 5(b)). Owing to the superposition theorem (Box 1, Eq. (2.3)), the superposed event sequence is itself a Poisson process with a rate $\Lambda$ that is equal to the sum of the individual rates, that is, $\Lambda =2\beta +\mu +\beta +\mu +2\beta +0=5\beta +2\mu$ .

Therefore, we can generate the time to the next event in the entire population using the inverse sampling method, which we introduced in Section 3.2, according to $\tau =-lnu/\Lambda$ , where $u$ is a uniform random variate in the interval $(0,1]$ . We now know the time to the next event but not which node (or more generally, which one individual Poisson process) is responsible for the next event. This is because the superposition lacks information about the individual constituent event sequences.

We thus need to determine which node produces the event. The mathematical properties of Poisson processes guarantee that the probability that a given node produces the event is proportional to its event rate (see Box 1, Eq. (2.4)). For example, in Fig. 6, where the event rates of the two nodes $v_1^{\prime }$ and $v_2^{\prime }$ are $2\beta$ and $\beta$ , respectively, each event in the superposed event sequence comes from $v_1^{\prime }$ and $v_2^{\prime }$ with probability ${\Pi }_1=2/3$ and ${\Pi }_2=1/3$ , respectively. Therefore, the probability that the first event is generated by $v_1^{\prime }$ is $2/3$ . The probability it is generated by $v_2^{\prime }$ is $1/3$ . This is natural because the sequence for $v_1^{\prime }$ has on average twice as many events as that for $v_2^{\prime }$ . In our example (see Fig. 5), the $v_i$ (with $i=1,\dots ,5$ ) that produces the next event is drawn with probability ${\Pi }_i$ , where ${\Pi }_1={\Pi }_5=2\beta /(5\beta +2\mu )$ , ${\Pi }_2={\Pi }_4=\mu /(5\beta +2\mu )$ , and ${\Pi }_3=\beta /(5\beta +2\mu )$ . We will explain computational methods for doing this later.

Figure 6 Superposition of Poisson processes and how to determine which component Poisson processes contribute to an event in the superposed event sequence. We consider $M=2$ Poisson processes, one with rate ${\lambda }_1=2\beta$ and the other with ${\lambda }_2=\beta$ . The superposed Poisson process has rate ${\lambda }_1+{\lambda }_2=3\beta$ . The next event in the superposed event sequence (shown in the dotted circle) belongs to process 1 with probability ${\Pi }_1={\lambda }_1/({\lambda }_1+{\lambda }_2)=2/3$ and process 2 with probability ${\Pi }_2={\lambda }_2/({\lambda }_1+{\lambda }_2)=1/3$ .

Assume that $v_4$ generates the next event and transitions from state I to state R. Then, the new event rate for each node is as shown in Fig. 2(b). We advance the clock by $\tau$ and go to the next step. Again, to determine the time to the following event, regardless of which node produces the event, we only need to consider the sum of the event rates of the six nodes, which is now given by $\Lambda =\beta +\mu +\beta +0+\beta +0=3\beta +\mu$ . Then, the time to the next event, which we again denote by $\tau$ , is given by $\tau =-lnu/(3\beta +\mu )$ , where $u$ is a new uniform random variate. We draw another random number to determine which of the four eligible nodes, $v_1$ , $v_2$ , $v_3$ , or $v_5$ , produces the event and changes its state. (Note that $v_4$ and $v_6$ are recovered, so they cannot undergo a further state change.) The four remaining nodes are each selected with probability ${\Pi }_1={\Pi }_3={\Pi }_5=\beta /(3\beta +\mu )$ and ${\Pi }_2=\mu /(3\beta +\mu )$ .

In this manner, we draw $\tau$ , determine which node produces the event, implement the state change, advance the clock by $\tau$ , and repeat. This is Gillespie’s direct method.

Let us take a look at another example, which is the SIR model in a population composed of $N$ individuals, in which everybody is directly connected to anybody else (called a well-mixed population; equivalent to a complete graph; see Fig. 7 for an example). In contrast to the general network case, each individual in a well-mixed population is indistinguishable from the others. In the well-mixed population, it is not the case that, for example, an individual has two neighbors while another has three neighbors (such as in Fig. 2); they all have $N-1$ neighbors. Recall that recovered individuals do not change their state again. Therefore, it suffices to consider $N_{\text{S}}+N_{\text{I}}$ event sequences and their superposition to determine the time to the next event and which individual will produce the next event.

Figure 7 Complete graph with $N=7$ nodes.

However, in the well-mixed population, we can make the procedure more efficient. Because everybody is alike, we do not need to keep track of the state of each individual.

In the case of a network, we generally need to distinguish between different susceptible individuals. For example, in Fig. 2(a), the susceptible nodes $v_1$ and $v_5$ are different because they have different sets of neighbors. Furthermore, $v_1$ has degree 2, while $v_5$ has degree 3. Therefore, if $v_1$ gets infected and changes its state, it is not equivalent to $v_5$ getting infected. So, we must keep track of the state of each individual in a general network. By contrast, in the well-mixed population, such a distinction is irrelevant. Everybody is adjacent to all the other $N-1$ nodes, and the number of infectious neighbors is the same for any susceptible individual, that is, it is equal to $N_{\text{I}}$ . It does not matter which particular node gets infected in the next event. The only thing that matters for the SIR model in the well-mixed population is the number of susceptible, infectious, and recovered individuals, which is $N_{\text{S}}$ , $N_{\text{I}}$ , and $N_{\text{ R}}$ , respectively.

It thus suffices to monitor these numbers. If an infection event happens, then $N_{\text{S}}$ decreases by one, and $N_{\text{I}}$ increases by one. If an infectious individual recovers, then $N_{\text{I}}$ decreases by one, and $N_{\text{R}}$ increases by one. Because the number of individuals is preserved over time, it holds true that

$N_{\text{S}}+N_{\text{I}}+N_{\text{R}}=N$(3.5)

at any time. Because anybody is connected to everybody else, any susceptible individual has $N_{\text{I}}$ infectious neighbors, so it gets infected at rate $\beta N_{\text{I}}$ . Because recovery occurs independently of a node’s neighbors, every infectious individual recovers at the same rate $\mu$ .

Based on this reasoning, we can aggregate the event sequences of the $N_{\text{S}}$ susceptible individuals into one even before considering which method we should use to simulate the SIR dynamics. (Therefore, this logic also works for the first reaction method and for rejection sampling.) Each event sequence corresponding to a single susceptible individual has the associated event rate $N_{\text{ I}}\beta$ . The superposed event sequence is a realization of a single Poisson process with rate $N_{\text{S}}\times N_{\text{I}}\beta =N_{\text{S}}{N}_{\text{I}}\beta$ . If an event from this Poisson process occurs, one arbitrary susceptible individual gets infected. Likewise, we do not need to differentiate between the $N_{\text{I}}$ infectious individuals. So, we superpose the $N_{\text{I}}$ event sequences, each of which has the associated event rate $\mu$ , into an event sequence, which is a realization of a Poisson process with rate $N_{\text{I}}\times \mu$ . If an event from this Poisson process occurs, then an arbitrary infectious individual recovers.

In summary, in a well-mixed population we only need to consider two coupled Poisson processes, one corresponding to contracting infection at a rate $\beta N_{\text{S}}{N}_{\text{I}}$ , and the other corresponding to recovery at a rate $\mu N_{\text{I}}$ . In a single step of the direct method, we first determine the waiting time to the next event, $\tau$ . We set $\tau =-lnu/(\beta N_{\text{S}}{N}_{\text{ I}}+\mu N_{\text{I}})$ , where $u$ is a uniformly random $(0,1]$ variate. Next, we determine which type of event happens, either the infection of a susceptible node (with probability ${\Pi }_{\text{S}\to \text{I}}$ ) or the recovery of an infectious node (with probability ${\Pi }_{\text{I}\to \text{R}}$ ). We obtain

$\begin{array}{cc}{\Pi }_{\text{S}\to \text{I}}& =\frac{\beta N_{\text{S}}{N}_{\text{I}}}{\beta N_{\text{S}}{N}_{\text{I}}+\mu N_{\text{I}}}=\frac{\beta N_{\text{S}}}{\beta N_{\text{S}}+\mu },\\ {\Pi }_{\text{I}\to \text{R}}& =\frac{\mu N_{\text{I}}}{\beta N_{\text{S}}{N}_{\text{I}}+\mu N_{\text{I}}}=\frac{\mu }{\beta N_{\text{S}}+\mu }.\end{array}$(3.6)

If an infection event occurs, we decrease $N_{\text{S}}$ by 1 and increase $N_{\text{I}}$ by 1. If a recovery event occurs, we decrease $N_{\text{I}}$ by 1 and increase $N_{\text{R}}$ by 1. In either case, we advance the clock by $\tau$ and go to the next step. We repeat the loop until $N_{\text{I}}$ hits 0.