Principled simulation of agent-based models in epidemiology

After more than a century of sustained work by mathematicians, biologists, epidemiologists, probabilists, and other experts, dynamic models have become a vital tool for understanding and describing epidemics and disease transmission systems. Such models fulfill a variety of crucial roles including data integration, estimation of disease burden, forecasting trends, counterfactual evaluation, and parameter estimation. These models often incorporate myriad details, from age and social structure to inform population mixing patterns, commuting and migration, and immunological dynamics, among others. This complexity can be daunting, so many researchers have turned to stochastic simulation using agent-based models. Developing agent-based models, however, can present formidable technical challenges. In particular, depending on how the model updates state, unwanted or even unnoticed approximations can be introduced into a simulation model. In this article, we present computational methods for approximating continuous time discrete event stochastic processes based on a discrete time step to speed up complicated simulations which also converges to the true process as the time step goes to zero. Our stochastic models is constructed via hazard functions, and only those hazards which are dependent on the state of other agents (such as infection) are approximated, whereas hazards governing dynamics internal to an agent (such as immune response) are simulated exactly. By partitioning hazards as being either dependent or internal, a generic algorithm can be presented which is applicable to many models of contagion processes, with natural areas of extension and optimization. Author summary Stochastic simulation of epidemics is crucial to a variety of tasks in public health, encompassing intervention evaluation, trend forecasting, and estimation of epidemic parameters, among others. In many situations, due to model complexity, time constraints, unavailability or unfamiliarity with existing software, or other reasons, agent-based models are used to simulate epidemic processes. However, many simulation algorithms are ad hoc, which may introduce unwanted or unnoticed approximations. We present a method to build approximate, agent-based models from mathematical descriptions of stochastic epidemic processes which will improve simulation speed and converge to exact simulation techniques in limiting cases. The simplicity and generality of our method should be widely applicable to various problems in mathematical epidemiology and its connection to other methods developed in chemical physics should inspire future work and elaboration.

Stochastic simulation of complex mathematical models is a vital tool for understanding 2 and describing disease transmission systems. While early efforts by probabilists Bailey 3 and Kendall [1,2], and the biochemist and physician-epidemiologist pair Kermack and 4 McKendrick [3,4] were highly mathematical in nature, they were at best, coarse 5 approximations of natural processes. In the century hence, disease transmission models 6 have incorporated myriad details including age structure, commuting, migration, and 7 within-host immunology to better represent our understanding of how these systems 8 function. Complex simulation models have been used to great effect in understanding 9 rapidly changing epidemic situations, such as the outbreak of a novel pathogen, or 10 events which require immediate response, such as the 2001 veterinary epidemic of hand 11 foot and mouth disease in the UK [5]. When addressing urgent public health crises, 12 complex models of epidemic processes can be invaluable tools, serving as platforms for 13 data integration [6], estimation of current burden [7], forecasting trends [8], evaluation 14 of intervention strategies and counterfactual scenarios [9], among other roles. 15 These models are often formulated as compartmenal models and simulated as 16 stochastic jump processes, where a jump is a change in an individual's epidemiological 17 state. Sample paths (trajectories) are piecewise constant, and change only due the to 18 occurrence of discrete events that change the epidemiological state of individuals 19 (jumps). When jumps are only allowed to occur after exponentially distributed intervals, 20 the process is a continuous-time Markov chain (CTMC). Because the modeler is free to 21 design the states, state transitions, and associated distribution of inter-jump intervals as 22 they see fit, this class of models can be directly specified from the results of survival 23 analysis [10], allowing close alignment to empirical data and standard survival models. 24 In general, jump processes are also valued for their deep connection to deterministic 25 models (e.g., compartmental models) in limiting cases which can aid model verification; 26 such results are well known for CTMC models [11,12] and exist for some non-Markovian 27 extensions [13]. Furthermore, these models benefit from over a half-century of rigorous 28 mathematical and algorithmic study [14], especially in chemical kinetics, physics, and 29 operations research communities, which has led to a plethora of publicly available 30 algorithms to sample trajectories from such processes, as well as techniques for 31 statistical inference and model fitting. 32 The complexity of these models, however, frustrates analytic approaches and can 33 even thwart straightforward application of classic simulation techniques, making quick 34 development and application challenging, especially in epidemic response situations. 35 Incorporating non-Markovian dynamics is still difficult in most simulation frameworks, 36 and because the majority of "industrial strength" simulation software is designed for 37 relatively straightforward chemical reaction networks [15], peculiarities of 38 epidemiological simulation such as highly nonlinear force of infection terms, age and 39 location-based mixing, immunological dynamics, and other elaborations make direct 40 utilization of these software difficult. While there exist some open-source software for 41 simulation of epidemiological dynamics, many frameworks are restricted to simple straightforward way to turn a whiteboard description of a complex system into useable 48 code, and are a viable alternative to other methods of representing a model. Many 49 ABMs are developed as bespoke programs for a specific analysis, but the technical 50 complexity of implementation means that a variety of subjective decisions may be made 51 when writing simulation code. Not all design choices will lead to simulation algorithms 52 that necessarily have a limiting interpretation as some stochastic jump process, valuable 53 for both model verification and model modification based on the formal rules of the 54 stochastic process. 55 Here we describe a method to construct approximate ABM representations of 56 stochastic models in which agents with arbitrarily complex internal dynamics interact 57 through discrete events in continuous time. Our method relies on specification of a 58 discrete time step of size ∆t, over which interactions between agents are approximated 59 (dependent events); dynamics within an agent (internal events) are still simulated 60 exactly in continuous time. A significant contribution in this work is presenting a 61 generic algorithm to simulate systems that are relevant to a wide class of 62 epidemiological models, via approximation of dependent events which can help speed up 63 even highly complex ABMs. We also demonstrate that our algorithm approaches the 64 true continuous time jump process as ∆t → 0. To verify our method and provide 65 numerical comparisons to exact stochastic simulation, we use our algorithm to simulate 66 a Markovian and non-Markovian SIR (Susceptible-Infectious-Recovered) model. We 67 conclude with a discussion of strengths and weaknesses of our approach, as well as 68 fruitful next steps to generalize our method. We hope that our method gives 69 mathematical epidemiologists considerable freedom in designing a model to fit their 70 needs, and that by approximating dependent events, even highly complex models can be 71 feasible to simulate. We also expect our method will be of interest to researchers in 72 ecology, demography, and the quantitative social sciences.

75
Formally, a jump process is a stochastic process, X whose trajectories (sample paths) 76 are piecewise constant functions of time over a countable set of states, S so that 77 X(t) ∈ S, t ≥ 0. Here we restrict ourselves to considering time-homogeneous processes, 78 so that only the dwell time in the current state is relevant for the process. In order to 79 sample trajectories from X, we must specify hazard functions λ j (τ, s) associated with 80 each event j ∈ {1, . . . , M }. The hazard functions tell us the conditional probability of j 81 occurring in the next infinitesimal time interval [τ, τ + dt), if the process has dwelled in 82 state s ∈ S for some time τ (in the case where there is no dependence on τ , the hazard 83 will be a constant value and X is a CTMC). When j occurs, it is allowed to change the 84 state X in some way. Given a set of M events, simulation consists of sampling when the 85 next event occurs, which event was it, updating state appropriately, recalculating 86 hazards that change, and repeating this process [18]. The state space can be defined 87 implicitly, and can be (countably) infinite, so long as only a finite number of events 88 have non-zero hazards at any time.

89
To define the process that the ABM will sample from, we let X(t) = (s 1 , ..., s n ) be a 90 vector where n is the number of agents being simulated, and each s h is the state of 91 person h. Expanding state space in this was to achieve an agent-based representation is 92 known as disaggregation, and may be used to study ABMs via the technical condition of 93 lumpability [19]. Then, any events which affect person h and whose hazard function 94 depends on more elements of X than only s h is a dependent event; events which affect 95 person h and whose hazard function is only allowed to depend on s h is an internal event. 96 Consider the recovery event in an SIR model. The recovery event for person h does not 97 need to know about the states of anyone else in order to calculate the hazard of recovery 98 and so is an internal event. However, the infection event for individual h does require 99 knowing the state of other agents in order to compute the hazard, and therefore is a 100 dependent event.

101
To draw approximate trajectories from the ABM requires a choice of time step ∆t, 102 over which hazards for each agent's dependent events may only use information from 103 other agents at the start of the step, ignoring changes which occur during the time step. 104 Put another way, agents only exchange information at the start of each time step. Then 105 over a time step beginning at time t, a susceptible agent i is subject to infection hazard 106 (force of infection, hereafter FOI) β(u)I(t) for u ∈ [t, t + ∆t), where β(u) is the effective 107 contact rate [20]. Note that the contact rate could be allowed to depend on time, and 108 that the approximation is in setting the number of infected individuals to a constant 109 value (the number at time t). Dependent events can be simulated via rejection sampling. 110 Internal events, by definition will not depend on any other individuals and therefore will 111 not be approximated. This type of simulation is reminiscent of tau-leaping 112 techniques [21], but whereas in standard tau-leaping all hazard functions are 113 approximated as constants over the time step, our method only approximates dependent 114 events, and internal events are simulated exactly. 115 We postpone a complete description of our method to section Simulation Algorithm, 116 and first demonstrate that the accept-reject sampling we use for sampling dependent 117 event times can draw correct times in a base case.

119
Consider the a single susceptible agent who is subject to a single event, infection. Let 120 τ S→I be a random variable giving the time at which this agent becomes infected, tmax 121 be the end of the current time step, [tmax − ∆t, tmax ), and λ be the FOI which is valid 122 over that time step.

123
In our model, each agent stores its current state s h , the time at which it entered that 124 state tnow , and the next scheduled event time tnext and state s h . For all internal 125 events, simulation is identical to classic discrete event simulation techniques. While the 126 agent's next scheduled event time is less than tmax the agent will update their state  In the case of infection, a dependent event, if we naively sampled the time of 131 infection as τ S→I ∼ tnow + Exp(λ), we would be ignoring future stochastic changes in 132 FOI, which could change at tmax . Crucially, from the perspective of this agent, the FOI 133 is a stochastic quantity, because it depends on the states of other agents. To develop a 134 reasonable approximation, one needs to implement a rejection algorithm to sample 135 τ S→I .

136
The sampling is described graphically in Fig 1 A. During this time step, λ is constant 137 (red region), and the agent samples a putative time to infection τ S→I ∼ tnow + Exp(λ). 138 Note that there is no restriction tnow be equal to the start of the time step, because it 139 could have been a randomly sampled quantity from previous events. If τ S→I falls within 140 the remaining time for which that hazard is valid (purple region of size tmax − tnow, 141 which is to say τ S→I < tmax ), then we accept the sample and infect the agent at time 142 τ S→I . The probability of the putative time being accepted is 1 − e −λ(tmax −tnow ) .

143
If however we reject the putative time ( τ S→I ≥ tmax ) then we set the agents next 144 event time to be tmax , and do not change the agent's state. Then, on the next time 145 step, the agent sets their current time to the start of that time step and resamples τ S→I 146 using the newly updated FOI λ . The probability of acceptance is 1 − e −λ ∆t . Exponential random variate.

165
If instead of constant, the FOI is a deterministic piecewise constant function such 166 that λ n is the value of the FOI between [n∆t, (n + 1)∆t) similar reasoning applies but 167 the number of trials needed no longer follows a Geometric distribution. The time of 168 acceptance within the acceptance interval however, still follows a truncated Exponential, 169 because it is conditioned on the infection event occurring in that time step of (piecewise) 170 constant hazard. In this case, sampling τ S→I is equivalent to sampling the first event 171 time of an nonhomogeneous Poisson process (NHPP) with intensity function λ(t).

172
From [22], one way to sample from such a process is by inversion of the distribution 173 of inter-event times, which has distribution function F (τ ; λ) = 1 − e − τ 0 λ(u) du . To 174 sample the first event time, one should draw a uniform random number u between [0, 1) 175 and solve so that τ is the first event time: If the cumulative 176 hazard function Λ(τ ) = τ 0 λ(u) du is especially easy to invert, then by the random time 177 change (RTC) theorem, we can also sample from the distribution of t by sampling v 178 from a unit rate Exponential distribution and then solving τ = Λ −1 (v) [23].

179
When FOI is piecewise constant, Λ(τ ) will be piecewise linear so inversion of the 180 December 15, 2020 5/20 cumulative hazard will be the most straightforward sampling method. To compare the 181 accept-reject algorithm with inversion sampling of NHPP first event times, we 182 discretized a continuous intensity functionλ(t) = 1 8 sin (t−6)2π 24 + 1 8 . The function has 183 a period of 24 hours, with a maximum amplitude of 0.25 (corresponding to 1/4 events 184 an hour) with a vertical shift so that it is non-negative and a phase shift of 6 hours such 185 that the maximum intensity occurs each day at noon and minimum intensity at In Fig 2 we show that the accept-reject sampler for first event times of NHPPs is 188 equivalent to exact inversion methods given in [22,23]. In both panels the red curve is 189 the density function of first event times calculated directly from f (τ ) =λ(τ )e −Λ(τ ) , and 190 we drew 10 6 samples using each algorithm to construct the histograms. We see that 191 both the accept-reject algorithm and integrated hazard sampling are sampling from the 192 correct density function. consider what happens as we let ∆t → 0. Because the model is a set of simple counting 198 processes in continuous time, meaning that only jumps of size 1 are allowed, the 199 probability of two events occurring simultaneously is zero [24]. So when ∆t is near zero, 200 we expect that on each time step either 0 or 1 event will occur, regardless of how many 201 agents are in the system. With infinitesimal time steps, after a single event occurs, the 202 FOI for all agents will be updated immediately and the algorithm becomes exact.

204
Our agent-based simulation algorithm for stochastic epidemic models where each agent 205 is subject to a single dependent event, infection, is given in pseudocode below.

206
Generalization to multiple dependent events (such as multiple strains or routes of 207

257
To illustrate use of our agent-based simulation method, we simulate a Markovian and In Eq (1) we present the Kolmogorov forward equations (KFE) for the Markovian SIR 274 model, following [25]. Because much recent research into methods for solving KFEs 275 originates in statistical physics, which, unlike probability theory or mathematical 276 epidemiology commonly uses step operators, we present the equations using step 277 operators and work out a more familiar form, as presented in [26]. A brief introduction 278 December 15, 2020 8/20 to using step operators to simplify writing KFEs for stochastic jump processes is 279 available in S1 Appendix. The term P(S, I, R, t) is the probability for the system to be 280 in state (S, I, R) at time t, such that S + I + R = N . We consider a simple mass action 281 force of infection term to simplify the mathematics, so that the deterministic R 0 = β γ N , 282 but the method is not restricted to simple mass action.  [27,28]. For an up to date review of methods 288 for solving transition probabilities for common epidemic models, except the method 289 of [28], see [29].  To assess the ability of the agent-based model to sample from the correct probability 306 distribution over future states when simulating trajectories, we initiated a simulation 307 with initial conditions S(0) = 60, I(0) = 10, R(0) = 0, γ = 1 /3.5, and R 0 = 2.5. We

Effect of Time
Step on Accuracy 320 We expect that as ∆t increases, the accuracy of the approximate ABM will deteriorate. 321 We therefore evaluated the distribution P(S, I, t = 5|S(0), I(0)) using the same probabilities, we used [31] to solve Eq (1).

329
For each time ∆t, we calculated absolute error as the sum of differences between the 330 Monte Carlo estimate of the transition matrix and the exact transition matrix.

339
An alternative to computing transition probabilities from Eq (1) for checking that our 340 ABM is sampling from the correct process is to compare the final epidemic size 341 distribution computed by Monte Carlo simulation to an exact analytic result. For SIR 342 models, a closed form final epidemic size distribution was developed in [32,33] and 343 recently reviewed in [34]. This closed form solution for the distribution of final epidemic 344 sizes is particularly valuable because it will be affected by the distribution of duration of 345 infectiousness, meaning it provides a complete check on the the ability of a sampling 346 algorithm to simulate the SIR model.
(2) Because the equations use the MGF of the infectious period distribution, Eq (2) can 354 compute final epidemic size distributions for both the Markovian and non-Markovian 355 SIR model. 356 We compared the final epidemic size distribution sampled from the exact simulation 357 algorithm (MNRM) and the ABM with ∆t = 0.01 to the exact closed form probabilities 358 calculated from Eq (2) in Fig 6. We used initial conditions of N = 50 and m = 1, and 359 sampled 10 4 epidemic sizes from each stochastic simulation. The effective contact rate 360 was calculated to give R 0 = 2.5, and γ = 1 /5. The exact probabilities are given by red 361 dots, and the empirical probabilities are given by the purple horizontal lines; pointwise 362 95% confidence intervals were computed for each empirical probability by Wilson's score 363 method and the coverage interval is given as the shaded rectangle around the empirical 364 probability [35]. In both cases the empirical distribution from stochastic simulation is

368
In Eq (3) we present the Kolmogorov forward equations for the non-Markovian 369 (semi-Markov) SIR model, where the infectious period τ ∼ F is a random variable that 370 has a density function f = d dτ F , which may differ from an Exponential distribution.

371
Because we assume that infection events still occur according to the points of a Poisson 372 process, the contribution to the KFE from infection is the same as Eq (1).

373
However the contribution from the recovery term is more complicated. In particular, 374 we need to introduce the two-time joint probability distribution 375 P(S, I, R, t; S t−τ , I t−τ , R t−τ , t − τ ), which is the joint probability of the system being in 376 December 15, 2020 12/20 state (S t−τ , I t−τ , R t−τ ) at a time t − τ , and in state (S, I, R) at a later time t. In 377 Eq (3), the recovery term must sum over all possible states that the process could have 378 been in τ units of time in the past, where τ ranges from [0, ∞), and is weighted by 379 f (τ )dτ , the probability of an infectious period of duration τ . This means when a S 380 particle becomes an I particle at time t − τ it immediately samples a time to recovery 381 according to f . Those recovery events which complete at time t have probability f (τ )dτ 382 of requiring that amount of time to do so.

383
In this sense, the recovery events in the system are still controlled by the points of 384 the Poisson process generating infection events, but represent delayed effects of that 385 event. An excellent description of stochastic systems with delayed effects is given in [36]; 386 another slightly different approach to these systems is described by [37]. and second terms represent infection events that will push the state into, and out of 389 (S, I, R), respectively. The third term represents infection events that fired some time τ 390 in the past and whose delayed recovery event will push the process into state (S, I, R) 391 at time t. Likewise the final term represents those infection events whose delayed 392 recovery events are just about to push the process out of that state.

393
To draw exact samples from the non-Markovian SIR model, we used the MNRM 394 from [30] for arbitrarily distributed delays. We modified Algorithm 7 such that the 395 delayed reaction channels store completion events as a priority queue implemented as a 396 binary heap, after drawing the delay, τ from the appropriate distribution.  Figure 7. We visualized the bivariate probability distribution over 409 P(S, I, t = 5|S(0), I(0)) using contours to represent curves of constant probability. We 410 observed a very good equivalence between the transition probabilities sampled via the 411 approximate ABM (dashed contours) and exact sampler (solid contours).

430
We have presented a practical approach to simulating ABMs with attractive properties. 431 One interesting consequence of using a time step is that integration with other models 432 becomes easier. If certain hazards depended upon, say a mosquito model with some 433 discrete time step, that information could be exchanged between models at each step.

434
As the model is formulated with continuous hazards, the size of a time step can vary 435 over a simulation run, potentially widening applicability of the method. Furthermore, 436 its generality means there are several interesting avenues for further development.   [38] for an example), it is not clear what is the 456 best way to deal with separation of time scale in the ABM. A naive approach would 457 involve slower interaction terms being updated only at some integer multiple of ∆t, but 458 further work should be undertaken to characterize the best way to choose such an 459 updating scheme suitable for multiple time-scales, with the caveat that appropriate 460 schemes would be highly problem dependent.

461
Our method is not restricted to pure jump stochastic processes. One generalization 462 of a pure jump process is to assume that the system can be described by coupling 463 continuous state variables to the discrete variables such that between jump times, the 464 continuous part evolves according to a differential equation. Such systems are known as 465 piecewise deterministic Markov processes (PDMP) when the distribution of inter-event 466 times is given by an nonhomogeneous Poisson process whose intensity may depend on 467 the continuous variables [39]. This model representation may be natural to simulate 468 multi-scale models, where within-host immune dynamics follow a differential equation 469 model, like the method proposed in [40] to model host-pathogen interaction. In the case 470 where the deterministic dynamics are fully internal to each agent the implementation 471 simply involves solving a differential equation for each agent to sample their next state 472 and time. More complex internal dynamics, such as allowing the continuous state to 473 follow a stochastic differential equation, can be implemented in the same way, as long as 474 the next jump time for each agent's discrete state is computable from the solution of the 475 process (for example, from a first hitting time). These models can be useful to integrate 476 results from medical survival analysis [41].

477
The computational methods for epidemiological problems have borrowed heavily 478 from other computational disciplines, but some of the unique features of epidemiological 479 systems would benefit innovation and algorithms to addressing the kinds of problems 480 that arise. In particular, for simulation of complex epidemiological models, agent-based 481 models are a useful alternative to other modeling techniques, and can be both flexible 482 and extensible. In this work we have presented a generic algorithm to sample 483 trajectories from agent based models, for agents which may experience dependent or 484 internal events, parameterized by hazard functions. Approximation of a subset of event 485 hazards, those that depend on multiple agents' state, speeds up simulation. A key 486 contribution of our ABM simulation technique is that it will converge to a 487 continuous-time stochastic jump process in the limit of small step sizes. Our intention is 488 not to argue that models without a limiting interpretation are wrong, but that in many 489 practical cases being able to connect the ABM to a well defined stochastic process is  [42,43]. We hope that the simplicity of the 495 method presented in this paper can help researchers respond more quickly to urgent 496 situations where stochastic models are required.

March 2020
Step operators are a convenient way to simplify the expression of Kolmogorov forward equations (master equations; henceforth KFE) for jump processes. We introduce them here in the context of continuous-time Markov chains, following the exposition of Toral and Colet 2014. Mathematically, the step operator E is a linear operator on a function f (n) where n ∈ Z. Because the probability flux on the right hand side of the differential KFE has contributions from all states with nonzero transition probabilities, and states are expressed as vectors of integers, the step operator E can be used to write KFEs. The main advantage of using step operators to write the KFE is that one no longer needs to keep track of what specific elements of the probability distribution the flux is flowing in and out of, one simply needs to keep note of how many "particles" of each type (corresponding to specific positions in the state vector) are being created and destroyed for each event that is allowed to change system state.
The step operator E is defined as below; for events which may cause the creation or destruction of more than one particle, note that E is defined for powers l. Most important for actually writing KFEs is the final relation, (E l − 1), which simplifies the inbound and outbound probability fluxes into a single term.
When writing KFEs with the step operator defined in Equation 1, note that to represent an event which creates l particles, one would use E −l − 1 [f (n)]. This is because the probability flux will include an incoming component from the state with l fewer particles, and the outbound flux from the current state.
Finally, for systems with more than one type of particle (the state space is a vector with more than one element), we can compose step operators to represent events which simultaneously change multiple particle counts simultaneously. Examples include predation in the stochastic Lotka-Volterra model, where predation simultaneously increases the count of predators and decreases the count of prey, or infection in simple epidemic models, which decrease the count of susceptibles and increases the count of infecteds. For an event which changes the count of type 1 and 2 particles by l 1 and l 2 , one would use Equation 2. To evaluate E l1 1 E l2 2 [f (n)], one applies the operators on f (n) one after the other.
To illustrate use of the step operator E, we demonstrate their use on a Lotka-Volterra model. In this model there are two types of particles, prey X and predators Y . There are three events, which can be written using notation from chemical kinetics as follows: These correspond to reproduction of prey, predation, and death of predators. At any point in time, the probability of the state n 1 prey (X) and n 2 predators (Y ) at time t is denoted as P (n 1 , n 2 ; t). Reproduction of prey has per-capita rate β. Predation occurs according to simple mass-action with reaction constant γ. Finally, death of predators occurs with per-capita rate µ. Now we write out the KFE and repeatedly apply the step operator until we get back to the familiar componentwise form.

March 2020
Another metric to check that our ABM is sampling from the correct process is to compare the final epidemic size distribution computed by Monte Carlo simulation of the ABM to an exact analytic result. We rely on the method presented in ) and well described in (Andersson and Britton 2012) to compute final size distributions for SIR epidemic models with general distributions of infectious periods.
Next pull out the final summation of the left hand side: Now multiply all terms through by Ψ F − (N −k)λ N k+m , noting that we can move the power in the denominator up to the numerator by adjusting the power by (k + m) − (i + m) = k − i: Finally, solve for p (N ) k , the desired probability: In order to turn Equation 6 into an algorithm, we note that when k = 0 the sum drops out and the equation simplifies to: Probabilities for k = 1, . . . , N can then be solved for recursively. It should be noted that for N much larger than 100, the algorithm given by Equation 6 suffers from numerical instability, and so should only be used for small populations. Because the equations are defined solely in terms of the MGF of the infectious period distribution, the probabilities from Equation 6 can be used for both the Markovian and non-Markovian SIR model.