STOCHASTIC SIMULATION ALGORITHM FOR EFFECTIVE SPREADING DYNAMICS ON TIME-EVOLVING ADAPTIVE

. Modelling and simulating of pathogen spreading has been proven crucial to inform containment strategies, as well as cost-eﬀectiveness calculations. Pathogen spreading is often modelled as a stochastic process that is driven by pathogen exposure on time-evolving contact networks. In adaptive networks , the spreading process depends not only on the dynamics of a contact network, but vice versa, infection dynamics may alter risk behavior and thus feed back onto contact dynamics, leading to emergent complex dynamics. However, numerically exact stochastic simulation of such processes via the Gillespie algorithm is currently computationally prohibitive. On the other hand, frequently used ‘parallel updating schemes’ may be computationally fast, but can lead to incorrect simulation results. To overcome this computational bottleneck, we propose SSATAN-X. The key idea of this algorithm is to only capture contact dynamics at time-points relevant to the spreading process. We demonstrate that the statistics of the contact-and spreading process are accurate, while achieving ∼ 100 fold speed-up over exact stochastic simulation. SSATAN-X’s performance increases further when contact dynamics are fast in relation to the spreading process, as applicable to most infectious diseases. We envision that SSATAN-X may extend the scope of analysis of pathogen spreading on adaptive networks . Moreover, it may serve to create benchmark data sets to validate novel numerical approaches for simulation, or for the data-driven analysis of the spreading dynamics on adaptive networks .


Introduction
Modelling and simulation of contagion spreading to forecast disease prevalence and to assess different public health interventions has a long history in mathematical biology [38]: Traditional approaches involve compartmental models, such as susceptible-infectious-recovered (SIR) or susceptible-infected-susceptible (SIS) [7].However, traditional compartmental models do not capture relevant spreading paths and may even provide incorrect predictions [37].Network models incorporate more realism, explicitly considering interactions (i.e.edges) between individuals and locations (i.e.nodes).Within this context, mostly static networks have been studied in the past [46,49].However, in static networks, the extent of linkage of nodes is determined only after integrating information derived from observing the contact (or spreading) process over a period of time.Several examples highlight that analysis of static networks lacks important temporal information about causal paths that underlie the spreading process, consequently yielding false conclusions for the control of the spreading process [32].In order to capture the temporal causality of the underlying system, different time-evolving network models have been introduced recently [33,39].Some of these approaches, in particular multi-agent-based models (ABM) may even consider spatio-temporal movement driven by undirected- [47] or directed dynamics [12].In these approaches, agents are typically deemed 'in contact' when they are in close physical proximity.Hence, the spatial dynamics are used to define the evolution of a temporal contact network.A particularly relevant extension to model spreading on contact networks is the inclusion of adaptive behaviour [54], which defines the class of adaptive networks.In this class of network models, the dynamical rules that shape the network structure change in response to the dynamics of the spreading process [27,28].Examples are the concurrency of sexual partnerships which is thought to be important for HIV spread [25,40,41] or measures of self-isolation (quarantine) for individuals diagnosed with SARS-CoV-2.This creates a feedback loop between the spreading dynamics on the network and the dynamics of the network itself, leading to emergent complex behaviour.Epidemic spreading has been extensively studied on these types of networks to understand the influence of social contact structure on disease prevalence [4,22,36,50].However, the dynamics of contagion spreading on adaptive networks are usually complex, and analytical results can only be obtained in special cases [21,50].This makes numerical studies based on stochastic simulations indispensable.
Several different approaches are available to date: The stochastic simulation algorithm (SSA) [23] allows for the exact numerical simulation of contagion spreading on adaptive networks, as used in [5,11,34].Readyto-use implementations are available in the network modelling packages Largenet2 [61] and EpiFire [31].A major drawback of the SSA is that it may generate computational overhead, when the contact dynamics are considerably faster than the contagion dynamics, i.e. if only a tiny fraction of contacts result in contagion transmission.While the per-contact transmission probability can be of the order 10 −1 for some airborne diseases like influenza [44], it is usually much lower in sexually transmitted diseases like HIV (order 10 −3 −10 −2 ) [6,52].The computational overhead becomes even more pronounced, when the outcome of intervention strategies, like face masks, vaccination, prophylaxis or 'treatment as prevention' are studied which further lower the per-contact contagion transmission probability [15][16][17][18].
Two types of algorithms are used in this context: (i) Inexact methods discretize the time and then perform parallel updates of the contact structure and the contagion spreading (akin to a tau-leaping procedure in systems biology [9,24]).Time discretization and synchronous updates are implemented in EpiModel [35], NepidemiX 1  and CovaSim [39].(ii) The method proposed by Vestergaard [57] borrows ideas from the 'integral method' originally proposed and applied in the Systems Biology field [1].It assumes deterministic contact dynamics, e.g. in Vestergaard [57], a temporal set of static contact configurations is used.If the contact dynamics would in fact be deterministic for the given time intervals, the method then estimates the exact time to the next spreading event.However, in Vestergaard's original implementation the underlying contact network is not adaptive any longer; i.e., the network configurations affect the spreading dynamics, but not vice versa.Therefore, this approach is limited to particular scenarios, in which human behaviour is not adapted as the pandemic unfolds, in contrast to what is currently observed for, e.g. the SARS-CoV-2 or monkey pox outbreak.
In this article, we develop an efficient and exact stochastic simulation method that allows to predict contagion spreading in the case when both the contact dynamics, as well as the spreading dynamics, are governed by inter-dependent stochastic processes, i.e. human behaviour is adapted as a function of the epidemic process.The algorithm is written in C++ and available through https://github.com/nmalysheva/SSATAN-X.

Model: spreading process on an adaptive contact network
In the following, we introduce the proposed simulation algorithm using an example of a simple spreading process that evolves on a time-dependent, adaptive contact network.The running example will involve adaptivity in the form of discontinuous re-wiring dynamics as a numerically hard problem that is meant to test the proposed algorithm (SSATAN-X).Application of the proposed algorithm to other types of networks (Fig. 1: static, temporal or numerically easier to treat adaptive networks) are exemplified in the Discussion.
Consider a population Ω = {S, I, D} that consists of susceptible-S, infected-I and diagnosed individuals D. Each of these individuals has a number of contacts that change over time.Infection (e.g.transmission of a virus) is possible from an infected to a susceptible, or from a diagnosed to a susceptible individual, only if the respective individuals are in contact.For an airborne pathogen, e.g.influenza or SARS-CoV2, a 'contact' denotes a period of time where individuals are in close proximity (e.g.within the same room), whereas for sexually transmitted infections it denotes a sexual relationship.
The contacts change over time.To describe the contact dynamics of the model, each individual is assigned rates of gaining and loosing contacts over time.
Moreover, in case the individual is diagnosed with the infection, he/she may become aware of his/her status, which affects the individual's behaviour.In the modelling example, this is where we include the (discontinuous) adaptive behaviour.In our model, the individual cuts all of his/her contacts, and the individual rate of establishing new contact drops to 30% of its pre-diagnosis value.Thus, this simple model describes adaptive dynamics, where the contact dynamics change depending on the epidemic state in a discontinous fashion, Figure 1C.
The presented model can be described as an undirected weighted temporal Contact Network G(t) = {E(t), V (t)}, where -V denotes the set of nodes {v 1 (t), v 2 (t), ..., v N (t)} ∈ Ω (the infected, diagnosed and susceptible individuals).Each node v i ∈ V has the following individual characteristics: (i) a rate at which a new contact is made λ + i and (ii) a rate of loosing an existing contact Parameter γ represents the transmission rate from j to k, which is greater zero if node j is infected or diagnosed and node k is susceptible and otherwise equal zero.
We consider a continuous-time Markov process that evolves the contact network G(t) in time.The possible events that can occur on the network can be classified into two main groups: 1. Contact dynamics R = {r jk }: Assembling of a new contact.In our example, for each pair of nodes (v j , v k ), j = k which is not connected by an edge, the rate of assembling an edge is defined by λ + jk = λ + j • λ + k i.e. product of the assembling rates of the two nodes.Disassembling of an existing contact.For each pair of nodes (v j , v k ), j = k that are connected by an edge, the rate of disassembling is defined as λ the product of the disassembling rates of the two nodes.

Epidemic dynamics A = {a }:
An infection emanating from an undiagnosed, infected individual In many applications, we are only interested in the epidemic dynamics, e.g. to quantify the effect of particular epidemic control measures, such as contact restrictions or other behavioral changes, on the infection incidence.Figure 1.Methods for simulating spreading processes on networks.We consider two processes Y |Z : E → E and Z|Y : V → V .The first process Y |Z affects the edges E of the network (topology, contact dynamics), while the second process Z|Y affects the state of the nodes V (epidemic dynamics).A-C.Different approaches for modelling spreading processes on networks.A. In static network approaches, the network topology E, which determines the linkage and thus the spreading process Z, remains fixed.B. In temporal networks, the network topology evolves according to some contact dynamics Y : E → E, which is independent of the epidemic dynamics Z|Y .C. Adaptive networks denote a class of co-evolving stochastic processes, where both the contact structure affects the spreading process Z|Y , and the spreading process alters the contact dynamics Y |Z.D. Exact stochastic simulation [23] of spreading processes on adaptive networks.In the stochastic simulation algorithm, every change to the network is explicitly modelled, including changes in the edges/topology, as well as changes to the state of the nodes (empty circles vs. filled red circles).E. Proposed method: The proposed method only considers the net effect of topological changes on the epidemic dynamics.An upper bound B L is computed such that B L ≥ j∈Z a j (V, E), and the network topology is bulk updated for the time step ∆(B L ) ∼ exp (B L ) using τ -leaping.A change to the state of the nodes (epidemic process) is conducted if r ≤ a 0,Z /B L ; r ∼ U(0, 1) and a 0,Z = j∈Z a j (V, E).

τ(Zt
However, the transmission dynamics directly depend on the dynamics of the contact process.Moreover, an important characteristic of any adaptive network is that not only the contact dynamics affect the spreading process, but also vice versa (Fig. 1C).I.e. in the example above, the rate λ + i and number of contacts of the individual depends on the actual state of node i.In our example, we will sample λ + i (t) (1/day) and decrease it by a factor in case of a diagnosis.Moreover, all existing contacts are lost upon diagnosis.However, any time-dependent function is in principle possible (see Discussion).However, we chose these (numerically hard to treat) discontinuous changes to put the proposed algorithm to the test.

Simulation of spreading dynamics with SSA
The presented model can now be evolved using the stochastic simulation algorithm (SSA), Figure 1D.The SSA would sample the time to the next event, which could either be an event that relates to the contact dynamics, or an event that relates to the epidemic dynamics as outlined in Algorithm 1.If the main interest is in assessing 12 Execute chosen reaction the epidemic dynamics, this approach would create considerable computational overhead, particularly if the contact dynamics are considerable faster than the epidemic dynamics (since only few contacts lead to infection, the contact dynamics are always faster than the epidemic dynamics).
In the following, we will present an algorithm that is able to precisely sample the time to the next infection event, while performing bulk updates on the contact dynamics.Notably, the bulk updates preserve all statistical properties of the contact dynamics that are relevant to the spreading process (epidemic dynamics).Hence, the algorithm allows to study the relevant dynamics at full resolution, while reducing computational burden.

Effective simulation of spreading dynamics with SSATAN-X
In the proposed algorithm, we are splitting the above described stochastic process into a process Y comprising all reactions R belonging to the contact dynamics and a process Z that is driven by reactions A (epidemic dynamics).We then leap-forward the contact process (Fig. 1E), until a reaction of the epidemic process happens (Algorithm 2).Algorithm 2 implements the core idea of SSATAN-X and is an adaptation of the EXTRANDE method proposed by Voliotis in the context of Systems Biology [58].
The original EXTRANDE algorithm was proposed for simulating "well-stirred" chemical systems, where each chemical particle ('agent') can interact with any other 'agent'.Our contributions involve the transfer of the algorithmic ideas to an evolving, sparse dynamical contact network that determines which agents are in contact with one another at any given time.Our contributions to enabling the application of EXTRANDE to aforementioned systems are fourfold: (i) We prove that EXTRANDE is an exact algorithm.(ii) We provide a method to calculate the (contact-structure dependent) upper bound for the epidemic dynamics, which is necessary for the exactness of EXTRANDE in simulating spreading processes on adaptive networks.(iii) We introduce a concept, that combines both EXTRANDE and bulk updating of the contact dynamics to enable computational efficient simulation of such systems.(iv) We adapt the algorithm to adhere to constraints imposed by actual contact networks (e.g. the number of contacts that can be formed and dissolve over time).
Algorithm 2: SSATAN-X envelope algorithm for simulating effective spreading dynamics on adaptive networks.
initialize: Network G(t), epidemic reaction functions a ∈ A, final time Define propensity upper bound B T L as proposed in Appendix S2: Bulk update the edges of G(t) for time interval [t, t + ∆t] using Algorithm 3 In the outline of the algorithm, line 3 defines the so-called Look-ahead time T L , for which an upper bound of the sum of propensity functions of the epidemic dynamics B T L for the network G(s) ; s ∈ [t ; t + T L ] is calculated.In Appendix S1 we proof that the algorithm is exact, if B T L ≥ a 0 (s) exists.According to [58] the look-ahead time horizon can safely be chosen as T L = T F − t.
Correspondingly, for T L , the upper bound B T L of the sum of reaction propensities of the epidemic dynamics is calculated in line 4.
Note that the parameter B T L defines an upper bound for the sum of all reactions propensities of A (epidemic dynamics) within the interval [t, t + T L ].It is therefore model-specific.

In our example, the set of reactions A consists of the transmission-(S → I), diagnosis-(D) and death reactions (∅). Therefore, we can define
B D,T L and B ∅,T L can be estimated straightforward as B D,T L = i δ i∈I and B ∅,T L = i∈I β i + j∈D β j .The upper bound B S→I,T L on the other hand depends on the amount of edges eligible for pathogen transmission, i.e. the number of S-I and S-D edges.In order to estimate the upper bound for the transmission rate, it is necessary to estimate the maximum possible number of edges/contacts that can lead to a transmission during the look-ahead time s ∈ [t ; t + T L ].As can be imagined, a large upper bound will result in many "thinning" reactions, which can make the algorithm inefficient.In Appendix S2, we provide a general procedure to calculate B T L and a tailored method to calculate a tight upper bound for the considered model.With respect to the upper bound B T L we can propose the time to the next epidemic reaction, ∆t (line 5).Notably, since B T L is an upper bound for the sum of the epidemic reactions a 0 (s) for network G(s) ; s ∈ [t ; t + T L ], the proposed time to the next reaction might be shorter than the actual time to the next reaction event.This apparent under-estimation is "corrected" by the "thinning" step (lines [18][19], such that the statistics of the process are being preserved (for details refer to [58] and for the proof to Appendix S1).When the proposed time step ∆t lies within the look-ahead horizon T L , time t is updated by ∆t (line 10).At this point, the actual propensity functions a need to be updated.In order to do so, the contact dynamics (edges) of the network G(t) are bulk-updated using Algorithm 3 (tau-leaping)(line 9) and the next reaction to update (or "thin" ) the state of the nodes (epidemic process) is chosen (lines [13][14][15][16][17][18][19].Here, the adaptivity step (line 17) refers to any drastic changes to the contact network, when an epidemic reactions occurs.In our example, a diagnosed individual breaks all his/her contacts.The bulk-update (line 9) allows us to avoid the unnecessary computation of each individual contact update by reactions R of the contact process.

Bulk update of contact network
The bulk-updating algorithm is based on the τ -leaping algorithm, which was originally proposed by Cao and Gillespie [9] and was later modified by Anderson [2] by incorporating post-leap checks to guarantee accuracy.
We use the bulk-update algorithm only for reactions R of the contact dynamics.A key adaptation is to consider all reaction of one type ξ ∈ {+, −} (addition or deletion of edges) simultaneously.In the network example, only two types of reactions change the contact structure of the network: R + = {r + jk } denotes the set of reactions that generate new edges {v j , v k }, j = k with rate λ For a bulk network update, there is a limitation of how many edge 'addition' and 'deletion' reactions can be executed: The number of edges available for deletion is limited by the number of edges existing in the network |E(t)| at time t.In total, the number of edges that exist in a network cannot exceed N (N − 1)/2 for a finite population of N individuals.Therefore, the number of edges available for addition is limited by N (N − 1)/2 − |E(t)|, which can be readily computed by considering a complement network.
We denote a contact graph (network) G = {E, V }, as well as its complement graph G = {E , V }.To accurately capture the contact dynamics, we store both these graphs.The contact graph defines its complement: i.e., a deletion of an edge in the contact graph is accompanied by the addition of the same edge in the complement graph.The use of these data structures allows to apply less complex search algorithms to correctly simulate the contact dynamics.Adding an edge to G will be complemented by deleting this edge from G and vice versa.Thus, we can describe the reactions of assembling and disassembling edges by the following first order reactions: Edge addition: If reaction r + jk fires with rate e jk • λ + j • λ + k , set e jk = 0 and e jk = 1.
Algorithm 3: τ -leaping for bulk-updating the contact network initialize: Network G(t), contact dynamics reaction functions r ∈ R, final time T F , 0 < p < p * < 1 and 0 < q < 1. q * > 1, < 1 for ξ do  τ -leaping with post-leap checks.Algorithm 3 describes the τ -leaping algorithm, modified for bulkupdating of the contact dynamics of the network.Besides the contact graph G(t), the contact dynamics reactions and the time horizon, the algorithm requires a few design parameters p, p * , q and that can, in principle, take many values.We used the optimal values derived in [2,9], e.g.we set p = 0.75, p * = 0.9, q = 0.98, = 0.03 and q * = 1.02.
Based on the initial state of the system, the time step of the leap τ is calculated (line 2) based on the following formulas: where |E|, |E | denote the number of edges in the contact network and its complement graph and = 0.03 is a design parameter.The expected net change to the network graph, as well as the complement graph is given by: where r +/− 0 denote the sum of reaction propensities for addition and deletion respectively at the current network state G(t).The variance of the net change is given by because both addition and deletion reactions are first-order reactions in the contact network and its complement graph.If the chosen τ is less than the value of 10 r 0 , we abandon the leaping algorithm and execute a fixed amount of steps using the SSA algorithm (usually taken to be 100 [24] lines 5-7 Find first y such as: Remove edge e jk corresponding to index y from the network and remove the corresponding element from ς − Find first y such as: Add edge e jk corresponding to index y to the network and remove the element from ς +

Implementation and availability
We implemented the algorithm in C++ using an object-oriented approach.The contact network was implemented with the LEMON Graph library2 .

Results
Unless otherwise stated, in this section we will generally use the following setup: We initialized the population with 200 nodes {S = 180, I = 20, D = 0} and an initial amount of 3000 edges.Each node j was assigned a parameter λ + j for establishing a new contact and a parameter λ − j for loosing an existing contact.These parameters were drawn from uniform distributions within the range (0.5, 2.5) and (0.4, 2.0) respectively.For the epidemic dynamics, the transmission rate along an S-I edge (contact) was set to γ = 0.04 and it was set to γ/2 along an S-D edge.The death rate of infected and diagnosed individuals was set to β = 0.08 and for susceptibles it was set to β = 0.The diagnosis rate was set to δ = 0.5.In case of a diagnosis event, where the infected individual is becoming aware of his/her status, the newly diagnosed individual cuts all contacts and the individual's rate of establishing new contact drops to 30% of the pre-diagnosis level.

SSATAN-X accurately captures contact network dynamics
In order to analyse the accuracy and performance of SSATAN-X, we first evaluated simulations with pure contact dynamics.Each edge is drawn as a bar for the duration of its existence.A. SSA captures every single event (deletion or addition of an edge) individually, resulting in very small time steps.B. The bulk-updating in SSATAN-X allows to create and/or delete several edges using larger time steps, while maintaining the correct statistics at these discrete time steps.
Whereas the exact stochastic simulation algorithm (Algorithm 1) conducts changes to the contact dynamics network one-by-one, as illustrated in Figure 2 (left), SSATAN-X uses bulk updates of the contact network, as described in Algorithm 3 and shown in Figure 2 (right).We remind that the main principle of SSATAN-X is that it accurately captures the statistics of the contact dynamics that is relevant to modelling the epidemic dynamics.Therefore, the statistics of the contact dynamics only need to be accurate at the time point when a bulk update is performed, but not in between bulk updates (Fig. 2, right).
To assess the accuracy of SSATAN-X in capturing the statistics of the contact dynamics, we first ran pure contact dynamics simulations, i.e. in which no epidemic transitions occur; a = 0 for all epidemic reactions .
Figure 3 shows the degree distribution (mean ± standard deviation) of the contact network at different time points using SSA vs. SSATAN-X respectively.The simulation starts with a deterministic initial network, as can be seen by the absence of the shaded areas in Figure 3A.During simulation, the degree distribution of the temporal network changes (Fig. 3B-D), until it reaches a stable distribution (Fig. 3E-F).During the entire simulation interval, both the mean degree distribution (solid blue vs. dashed red lines), as well as the standard deviations (blue and red shading) remain visually indistinguishable between the two algorithms.
To further quantify the differences between the evolving contact network over time, we computed the statistics of the Kolmogorov-Smirnov test, as shown in Figure 4A-B.As can be seen in Figure 4A, the probability that the degree distributions generated with the SSA vs. SSANTAN-X are identical is much bigger than 0.95 (Pvalue) for all time points.The test statistic of the Kolmogorov-Smirnov test (distance between the respective ECDFs) is shown in Figure 4B.The test statistic remained below the value of 0.0015 at all time points, further highlighting that differences in the degree distribution are insignificant.
We hence conclude that the bulk-updating (Algorithm 3) of the contact network dynamics, as implemented in SSATAN-X, accurately captures the statistics of the time-evolving contact network at the discrete time steps.
Next, we evaluate whether the SSATAN-X also accurately captures epidemic dynamics that unfold on an adaptive contact network.The dashed red and solid blue lines represent mean degree distribution over 10 3 simulations using the SSA and SSATAN-X respectively.Shaded areas represent the corresponding standard deviations for the 10 3 simulations.

SSATAN-X accurately captures epidemic dynamics
In Appendix S1, we prove that Algorithm 2 is exact for simulating the epidemic dynamics if an upper bound B T L for the sum of propensities of the epidemic dynamics a (s) s ∈ [t, t + T L ] exists.In this section we present the results of the SSATAN-X algorithm applied to the adaptive susceptible-infected-diagnosed (SID) network model described in the Methods section.
Foremost, we were interested to see whether the statistics of the times to the next epidemic event were correct, as proposed in Figure 1D-E.In Figure 5, we show a histogram of the time-steps to the next epidemic event (infection, diagnosis, death) for the SSA (orange bars), as well as SSATAN-X (blue bars).As can be seen, both algorithms yield indistinguishable statistics regarding the epidemic inter-event times.This reaffirms the correctness of the proposed algorithm in being able to modelling the net effect of contact dynamics on the spreading process.Figure 6 shows the corresponding simulation results for the epidemic dynamics.The mean (± standard deviation) number of susceptible, infected and diagnosed are shown in Figure 6A-C respectively.Both the mean and the standard deviation are indistinguishable when comparing simulations with the SSA (red line and red shading) vs. SSATAN-X (blue line and blue shading) over the entire duration of the simulation.We used the Kolmogorov-Smirnov test to compare the degree distributions obtained by SSA and SSATAN-X (represented on Fig. 3). A. The P -value denotes the probability that the two degree distributions derived from SSA-and SSATAN-X simulations respectively are identical.B. Statistics of the Kolmogorov-Smirnov test (KST), representing the distance between the two empirical cumulative distribution functions.Therefore, we can further confirm the correctness of SSATAN-X in simulating epidemic dynamics on an adaptive contact network.

Impact of contact network adaptivity on the epidemic dynamics
The analysed S-I-D model (Methods section), belongs to the class of adaptive network models, because not only do the contact dynamics affect the epidemic dynamics, but vice versa, the epidemic dynamics also feed back onto the contact dynamics (compare Fig. 1C).Specifically, in the utilized model, the epidemic event of diagnosis changes both the contact network itself (cutting all edges of a newly diagnosed individual), but also its further dynamics (the rate of establishing new contacts drops to 30%).
To investigate the effect of adaptivity on the epidemic, as well as the contact dynamics, we consider two representative examples of the contact network model.For the first of the examples we set the diagnosis rate to δ = 0, in other words, in this example we have a non-adaptive temporal network model (described in Fig. 1B).For the adaptive model (second example), we set the diagnosis rate to δ = 0.5.In both examples, we set the transmission rate γ to 0.004.Adaptivity in the second example is implemented in the following way: The diagnosis event of an infected individual causes a complete rewiring of the individuals' contact network.The individual looses all contacts upon diagnosis and subsequently has a decreased rate of making new contacts, i.e. λ + j,diag = 0.3 • λ + j,undiag , where the subscripts 'diag' and 'undiag' refer to the individual being diagnosed and undiagnosed with the infection respectively.The transmission rate γ remains unchanged.Also, for both examples we set the death rate to β = 0.This allows to investigate the impact of adaptivity on the epidemic and contact dynamics, without the impact of other factors.
Figure 7A shows the number of edges in the network that are eligible for transmission in the non-adaptive temporal network (δ = 0; blue error bars), vs. the adaptive network (δ = 0.5; orange error bars).Annotations show the percent of infected individuals that are diagnosed (e.g.D/(D+I)).Figure 7B compares infection dynamics of both, the non-adaptive and adaptive model (blue vs. orange error bars).As can be seen, the diagnosis event alters the epidemic dynamics in the sense that less new infections occur in the adaptive network.In particular, since we set up the simulations in such way that the transmission rates emanating from both infected and diagnosed were identical, the decrease in infection dynamics can be solely attributed to alterations in the contact dynamics in the adaptive networks (Fig. 7A), which then feed back onto the epidemic dynamics.i.e., the number of edges, eligible for the transmission of the virus decreases with an increasing percentage of infected individuals being diagnosed (Fig. 7A).The consequence is fewer opportunities for disease spreading and hence the number of new infections (orange error bar in Fig. 7B) decreases, too.

SSATAN-X speeds up computation time
Next, we analyzed the run time of SSATAN-X in comparison to SSA for different populations sizes (number of nodes).To this end, we took the model settings and parameters as described in the beginning of the Results section.
Simulations were started with fixed number of infected individuals {I = 20}.The simulation results presented in Figure 8A show that the run times increase with increasing population size.However, for the utilized parameter setting, the runtime of SSATAN-X is about 100-fold less compared to the run time of SSA.
Next, we wanted to assess the runtime performance as a function of the utilized parameters.To that end, we chose homogeneous contact dynamics rates λ + i , λ − i across the population (all individuals having the same rates).This allows to assess how the run time changes as a function of the probability that transmission occurs  on a contact edge before it dissolves.This probability is given by γ γ + (λ − ) 2 , where γ is the transmission rate on a contact edge and (λ − ) 2 is the rate of loosing an existing edge.In the analysis, we increased λ + and λ − proportionally (the average number of contacts stays the same, but their lifetime shortens).
The run time of both SSA and SSATAN-X for performing 10 3 simulations for each parameter set is shown in Figure 8B.As can be seen, when the relation γ γ + (λ − ) 2 is close to 1, i.e. transmission is always bound to happen before an edge dissolves, SSA may be slightly faster than SSATAN-X.The reason is that the timeinterval between two epidemic events becomes so small that the bulk update of the contact dynamics reduces to the SSA.Noteworthy, such models would actually not require to model separate contact dynamics altogether, making this example a bit artificial (more on the utility of different modelling approaches in the Discussion section).
When the relation γ γ + (λ − ) 2 decreases, i.e. the contact dynamics being much faster than the epidemic dynamics, the run time of SSATAN-X becomes superior over SSA.Therefore, with low transmission rates and fast contact dynamics, SSATAN-X avoids any computational overheads from the contact dynamics and only takes their net effect on the epidemiological dynamics into account by bulk updating.On the contrary, SSA requires to execute each reaction (contact and epidemic) individually.
In summary, SSATAN-X is computationally more efficient than SSA.The run time of both algorithms scales with the population size and the "boost-factor" of SSATAN-X is preserved across different population sizes.However, the computational "boost-factor' of SSATAN-X increases, when the contact dynamics are fast in comparison to the epidemic dynamics.

Comparison with parallel updating methods
Many softwares for simulating spreading dynamics on networks have recently been introduced as a response to the ongoing SARS-CoV-2 pandemic, such as OpenABM-Covid19 [30] or CovaSim [39].Some of these softwares   2 , with γ being the transmission rate and λ − being the rate of edge disassembling, can be interpreted as the probability of transmitting an infection before the contact is dissolving.allow to consider adaptive networks (compare Fig. 1).For simulation, these softwares frequently use parallel updating schemes, where a constant time-step is pre-selected by a user or hard-coded in the software.In each time step, updates to both the contact network and epidemic state are performed.For example, the popular CovaSim software [39] uses a fixed time step of one day to perform parallel updates.
To compare SSATAN-X with these inexact methods, we implemented a parallel updating scheme (Appendix S3) and tested the accuracy of parallel updating schemes for different time steps on the previously introduced model.As the smallest time step for parallel updating, we used dt = 0.0001, i.e. smaller than the average timestep of SSATAN-X (dt = 0.062).i.e.: This should, on average, lead to less than a single epidemiological event and thus should lead to fairly accurate results with parallel updating.
As can be seen in Figure 9A-C, SSA (exact simulation) and SSATAN-X produce almost indistinguishable results with regards to the average number of infected, diagnosed and susceptible individuals (red and blue lines).Corresponding statistical ranges (mean ± one standard deviation) of the predictions are shown in Figure 9D-F.For the smallest considered time steps (dt = 0.0001 and dt = 0.1) parallel updating is relatively accurate, but nonetheless slightly overestimates the average number of infected (and diagnosed infected) individuals.This error propagates and it could be envisioned to worsen significantly over very long simulation horizons.The reason for this behaviour lies in the discontinuous contact dynamics that we considered: i.e. upon diagnosis, all edges are lost.In parallel updating these edges will remain lost until the next time step, whereas in reality (as modelled by SSA, SSATAN-X) some new (transmission eligible) edges may have arisen.This shows that SSATAN-X can even handle numerically very challenging dynamics and can consider events as they occur (for example: a person taking a diagnostic test just before making the decision to enter a large social event).
For larger time steps dt = 0.2, dt = 0.5, which imply that several epidemiological events are executed in parallel (approx.2-3 for dt = 0.2), the results from parallel updating become wrong altogether.This can be highly problematic for parallel updating schemes, which fix the time step throughout simulation: i.e. while the chosen time step in parallel updating schemes may be sufficient at the start of an outbreak simulation (few transmission occur), after some time, when the epidemic expanded, the same time step may imply several epidemic events, and this could lead to incorrect predictions.Thus, there is no guarantee for an arbitrary chosen time step, that simulations using parallel updating are correct for epidemics that unfold on adaptive networks.
For a more close inspection, we also computed the Kolmogorov-Smirnov Test (KST) statistics for the number of infected, diagnosed and susceptible individuals in Figure 10, compared to SSA.Note that this is assessment is different from Figure 4, which only assessed the contact dynamics (degree distribution).Here, we compare the statistical difference with regards to the epidemic states (susceptible, infected, diagnosed).As can be seen, the value of the KST statistics can be markedly higher for the parallel updating scheme than for SSATAN-X.For parallel updating with the smallest step sizes (dt = 0.0001), the KST statistic is of similar magnitude compared to SSATAN-X.For dt = 0.1, it is for the majority of depicted time points (0.5, 1, ..., 5) larger than the KST statistic of SSATAN-X (8/10, 9/10 and 5/10 for infected, diagnosed and susceptibles).For dt = 0.2 and dt = 0.5 the derived distribution is clearly different from SSA, hence simulation results become inaccurate.
For the utilized model and the smallest time step dt = 0.0001 parallel updates are about 10 times slower than SSATAN-X, which required about 1300 milliseconds for 1000 simulations on an Quad-Core Intel i7 laptop with 2,5 GHz.For larger time steps (dt = 0.1, 0.2 and 0.5), computational run times were half of that of SSATAN-X, at the cost of accuracy.

Discussion
Spreading of infectious disease can occur when individuals encounter in contacts that are relevant for the route of transmission of a pathogen.These contacts can be sexual contacts for sexually transmitted disease, or spatial proximity for airborne infections.In order to more realistically model spreading of infectious diseases, different network models have been proposed in the past [28].There are important circumstances in which the dynamic process of the contact network and the epidemic spreading are worth considering.Roughly speaking, this is when the two dynamic processes are on similar timescales.For example, if the contact dynamics are much slower than the epidemic dynamics, such that transmission is bound to happen whenever there is a  contact, it is sufficient to only consider the contact dynamics.Also, if the network changes so slowly that it is likely to remain the same throughout an outbreak, then a even static network approximation would be appropriate [19].However, if the network changes over a similar timescale to the node-to-node progression rate of the infection, then the changing network can change the course of the contagion, and it is important to include those changes in the analysis.A particular case, that we address with SSATAN-X is the case where the temporal processes of the contact network and epidemic interact not only in timescale, but directly: The infection events themselves change the network structure, which then impacts further disease spread.These more complex models, where a contact network structure changes dynamically in response to the dynamics of the spreading process, are often termed adaptive networks [27].These types of models may therefore allow to realistically consider behavioural changes that occur when individuals become aware of their infection and self-isolate [48], or experience debilitating symptoms.
In this article, we introduce SSATAN-X, an efficient and exact algorithm to simulate effective spreading dynamics on adaptive networks.
SSATAN-X takes advantage of the fact that for any model, where it is worth considering both contact dynamics and epidemic dynamics, the former must be faster than the latter (i.e.not every contact leads to transmission).Instead of performing detailed simulations of the underlying contact dynamics, the core idea of SSATAN-X is to only consider the effect that the underlying contact dynamics have on the epidemic dynamics.For this purpose, SSATAN-X uses bulk updates of the contact structure between epidemic events, such as disease transmission, diagnosis or death (Fig. 2).
In Figure 3, we show that the proposed bulk-updating conserves the statistical properties of the underlying contact network at times that are relevant to the epidemic dynamics.Because the statistics of the underlying contact network at times that are relevant to the epidemic dynamics are preserved, the statistics related to the firing of the next epidemic event (Fig. 5) are preserved, too.Consequently, the resulting epidemic dynamics are statistically correct, Figures 6 and 7.
However, the bulk updating avoids computational overheads from updating contact edges one-by-one.Therefore, depending on the choice of parameters, SSATAN-X can speed up computation time quite significantly in comparison to exact methods like the stochastic simulation algorithm (SSA) without any loss in simulation accuracy, Figure 8.The speed up factor critically depends on the relation of the contacts dynamics to the transmission dynamics.For example, if contacts generally do not result in transmission, computational speed up with SSATAN-X will be considerable.On the other hand, if transmission will always occur on a contact edge before it dissolves, the computational speed up will only be moderate or even absent.In Figure 8B, we expressed this dependency as γ/(γ + (λ − ) 2 ), where (λ − ) 2 denotes the rate of loosing an edge and γ is the transmission rate along a contact edge.
Note that the transmission rate γ denotes a lumped parameter, that can be further decomposed to allow for more sophisticated modelling, realistic parameterization and model extensions: For example, if exposures on a contact edge occur with rate r e and if the exposures are of the same type, then the transmission rate is just a scaled variant of the exposure rate, i.e. γ = r e • p tr|e , where p tr|e is the probability that a transmission occurs for a single exposure event e. Notably, for many infectious diseases and types of exposures p tr|e is readily known [6,44,52] and the effects of pharmaceutic [16,18] or non-pharmaceutic interventions [55,56] on p tr|e can be quantified, allowing for model extensions that can link within-host pathogen dynamics to the spreading dynamics at the population level.
Taken together, the lumped transmission rate γ is a function of the transmission probability per exposure event p tr|e , as well as the frequency of exposures r e when individuals are in contact.For many infectious diseases and types of contacts it may be argued that the ratio γ/(γ + (λ − ) 2 ) is much smaller than one; i.e., most transient contacts may not result in pathogen transmission.Behavioral changes, such as distancing, isolation or quarantine, may further reduce the transmission rate by reducing r e , whereas pharmaceutic (vaccination, prophylaxis, treatment) [10,26,51] or non-pharmaceutic interventions (e.g.mask wearing for airborne-, and condom usage for sexually transmitted infections) [42,59] may lower p tr|e in a static or time-dependent manner.
In our examples, we used a very simple model to demonstrate that the proposed algorithm works (proofof-principle).In the utilized model and in the algorithm to simulate spreading on this model we differentiate between (i) epidemic dynamics: For example, reactions that change the population (birth, death) or the state of the considered population (susceptible, infected, diagnosed) and (ii) contact dynamics: those reactions that affect the contact network of individuals in the populations.Note that in principle, the algorithmic idea applies to any model that can be set up this way, no matter which exact reactions belong to each of these kinds.With regards to the presented model, two adaptions are needed: the computation of the upper bound B T L and the computation of a (t) in line 10 of Algorithm 2.
In the examples used throughout, we considered reaction rates (for R and A) to be constants.However, it is also possible to consider time-dependent rates: For example, the transmission rate γ may be a function of time, due to some time-dependent pharmacokinetic effect of a prophylaxis [15,17,60], or depending on how long an individual is being infected or contagious [55].These time-dependent effects may be incorporated when determining the upper bound B in Algorithm 2 (Appendix S2).In fact, in Supplementary Figure S1, we extended the utilized model, by including a time-dependent force of transmission γ(t).Here, we used a 'sine' function, e.g. to represent seasonality, but any time-dependent function is possible.As noted above, this requires two adaptations: when computing B T L , the force of infection 'γ' is replaced by max s γ(s) for some s ∈ [t, t + T L ].After sampling the time step ∆t ∼ Exp(B T L ) and progressing to t ← t + ∆t, the actual γ(t) needs to be calculated, and based on it a (t).Note that the aforementioned example can readily reflect dynamic (time-dependent) infection control that affects the transmission rate, such as 'mask-wearing' for airborne infections, or vaccination and prophylaxis.If some time-dependent control of contacts is applied, this can be implemented two ways: In case of abrupt changes, i.e. contact restriction being on/off at a particular time, the implementation would be straight forward: The look-ahead time horizon T L is chosen such that it coincides with the change in the contact rates.At the rejection step (line 6 in Algorithm 2) the rates of establishing new contacts are changed accordingly, e.g. from the rates before to after contact restrictions.If the rate of establishing new contacts λ + (t) was an arbitrary time-continuous deterministic function (e.g. a 'sine' function or the like), the computation of B T L needed to be adapted.This can be done, by using the time dependent λ + (t) in equations ( 1) and (2) and equation (7) in the Appendix S2 "Determination of Upper Bound".The equations ( 1)-( 2) and equation ( 7) could be solved numerically using standard ODE solvers to calculate all variables in equation ( 15) Appendix S2 "Determination of Upper Bound".Note that if continuous-time functions were considered for either contact or epidemic reactions, the SSA algorithm (or parallel updating schemes) would be numerically inexact, since it is designed for simulating homogeneous Markov Processes.
Lastly, high-resolution networks representing e.g. an agent's mobility in the physical space can also be simulated using the proposed algorithm.Spatial movement and proximity can be modeled by having individual contact rates for pairs of nodes j,k.For example: In the manuscript, we considered reaction r + jk fires with rate e jk • λ + j • λ + k .i.e., we did not regard spatial proximity or movement and any two nodes can, in principle, come into contact.We did this, because we focus on presenting the basic algorithm, rather than complex modelling.However, if instead r + jk fires with rate e jk • λ + jk , where λ + jk considers whether the two individuals j,k meet on their way to work or work in the same premises, spatial proximity can be modeled.Note that the parameters for the contact dynamics may, in addition, be time-dependent, as outlined above to reflect complex scenarios, e.g.: 'individuals j and k are likely to meet between 9am and 10am on their way to work and not thereafter or before'.
Interestingly, if these spatial contact dynamics are NOT adaptive they may actually be pre-computed and plugged into Algorithm 2 to directly compute B T L : e.g., B T L (S → I) = γ • i∈I j∈S T L 0 δ ij (s) ds, where δ ij (s) takes the value 1 if an an edge exists at time s.
In our paper, we use the Anderson tau-leap algorithm (Algorithm 3) [1] to compute the bulk updates of the contact dynamics between two epidemiological events.It is also possible to use other tau-leap based approaches, as summarized in [45], however, we used the Anderson tau-leap algorithm, because it does guarantee accuracy through step size adaptation.
While τ −leaping may be considered as a stochastic simulation algorithm that has some analogies to explicit (first order) Euler methods for solving systems of ordinary differential equations (ODE) [8], another possibility for the bulk updating of the contact dynamics is to adapt higher-order schemes with adaptive step sizes that have been proposed in the context of ODEs [53].For example, the Runge-Kutta-Fehlberg method [20] could be used to adapt step sizes for the bulk updating.A clear advantage of higher-order schemes in the context of bulk-updating the contact network is the fact that larger step sized may be chosen that preserve accuracy [3].In the future, we will further improve SSATAN-X by incorporating these higher-order approximation schemes.
In our paper we consider finite populations.In finite populations, the maximum possible amount of contacts of the each individual i has an upper theoretical limit of N − 1 and a complete graph would have N • (N − 1)/2 edges.This gives us a straightforward way to estimate the total number of edges available for deletion or addition at any time point t: The number of edges available for deletion equals the number of edges present E(t), whereas the number of edges available for addition is given by N • (N − 1)/2 − E(t) or in other words, the number of edges in a complement graph.Note that for infinite populations the number of edges that can be deleted is still given by E(t), whereas, for adding edges one may simply abandon the 'complement graph' approach, making the algorithm more efficient.Some authors have also proposed more restrictive models for the number of contacts [40].Some of these concepts may be inspired by the popular Dunbar number [13,14], which is rooted in the presumption that individuals may not maintain stable social relationships with more than 150 other individuals.Notably, Dunbar's presumption has not been without critique [43].
However, contact networks may exist, in which the maximal number of contacts of an individual is less than the number of available individuals in the network.For example, Schmid and Kretzschmar [40] describe an individual-based model, in which each individual has an upper-limit capacity regarding the number of concurrent sexual partnerships, with pair formation and separation simulated as a dynamic processes, where each individual has his/her own probabilities of contact formation and separation drawn from an appropriate distribution.To model these kinds of restrictive contact dynamics with SSATAN-X, the number of edges available for deletion remains unaffected, while the estimation of the number of edges available for addition becomes more involved: Depending on the order of the addition and remaining contact capacities, there may be several outcomes with different amounts of edges possible to be added.For example, in a network where all nodes but one have exhausted their capacities, no more edges can be added without violating the contact constraints.Therefore, it becomes necessary to estimate a "worst scenario" case, i.e. a contact bulk update that resulted in the least

2 −
occurs with rate γ greater zero if node j and k are connected.An infection emanating from a diagnosed individual D i + S k γ/−→ D i + I N I +1 occurs with rate γ/2 greater zero if node j and k are connected.An infected individual may be diagnosed with the infection I j δ −−→ D N D +1 , for all j = 1...N I .Individuals may die:

Algorithm 1 : 4 Compute the sum of the reactions propensities 5 r 0 = r jk 6 a 0 = a 7 ∆t ∼ Exp (a 0 + r 0 ) 8 9 r jk r 0 + a 0 , 10 or an epidemic dynamics with probability: 11 a r 0
SSA initialize: Network G(t), reaction rate functions r ∈ R, a ∈ A, final time T F 1 t = 0 2 while t < T F do 3 Compute reaction propensities r (•,•) and a (•) based on the current contact network G(t)Choose the next reaction, e.g. a contact dynamics reaction with probability: + a 0 .

end 20 end 21 if
leap conditions equation (2.4) holds then for ξ do Delete all rows from S ξ less than or equal to row ξ Add new first row [r ξ 0 would have failed the leap condition for = 0.75 then τ = τ p * else Set τ = τ q if t < 1 and otherwise τ = τ q * end Bulk-update network by executing M ξ reactions using Algorithm 4 Recalculate r ξ (k,l) , r ξ 0 and r 0 based on updated network state 36 else for ξ do Add row [r ξ 0 τ + Q ξ , C ξ + M ξ ] between row ξ and row ξ + 1 end τ = τ p 41 end 42 end 43 end Edge elimination: If reaction r − jk fires with rate e jk • λ − j • λ − k , set e jk = 0 and e jk = 1, where e jk ∈ E ∈ G and e jk ∈ E ∈ G .

8
Calculate tmp = r + jk + value of last element in ς + 9 append tmp to ς + 10 else if addition then 11

14
Calculate tmp = r − jk + value of last element in ς − 15 append tmp to ς − 16 end

Figure 2 .
Figure2.Edge activity plot.Edges of the network depicted in the order of appearance.Each edge is drawn as a bar for the duration of its existence.A. SSA captures every single event (deletion or addition of an edge) individually, resulting in very small time steps.B. The bulk-updating in SSATAN-X allows to create and/or delete several edges using larger time steps, while maintaining the correct statistics at these discrete time steps.

AFigure 3 .
Figure 3. Degree distributions during simulation.The graphic depicts the mean degree distribution of the network during distinct time points within the simulation time (t ∈ [0, 5]).The dashed red and solid blue lines represent mean degree distribution over 10 3 simulations using the SSA and SSATAN-X respectively.Shaded areas represent the corresponding standard deviations for the 10 3 simulations.

Figure 4 .
Figure 4. Statistical assessment of differences in the simulated contact network.We used the Kolmogorov-Smirnov test to compare the degree distributions obtained by SSA and SSATAN-X (represented on Fig.3). A. The P -value denotes the probability that the two degree distributions derived from SSA-and SSATAN-X simulations respectively are identical.B. Statistics of the Kolmogorov-Smirnov test (KST), representing the distance between the two empirical cumulative distribution functions.

Figure 5 .
Figure 5. Time-step to the next epidemic event.The histogram depicts the frequency of the time-step sizes to the next epidemic event (infection, diagnosis or death), produced from 10 3 simulations for each algorithm and normalized over the number of samples (overall 79375 epidemic time steps for SSATAN-X (blue bars) and 79025 for SSA (orange bars)).

Figure 6 .
Figure 6.Simulated infection dynamics.The graphic depicts the change in the number of infected (panel A.), diagnosed (panel B.) and susceptible (panel C.) individuals over simulation time.Dashed red and solid blue lines describes the sample means over 10 3 simulations using SSA and SSATAN-X respectively.Shaded areas represent the corresponding areas encompassed by the mean ± standard deviation.

Figure 7 .
Figure 7. Impact of diagnosis on the degree distribution and amount of new infections.Blue error bars portray the results from a non-adaptive model, where we set the diagnosis rate δ = 0. Orange bars depict predictions from an adaptive model, where the diagnosis rate δ = 0.5.In the adaptive model, a diagnosis of individual j, leads to a loss of all contacts (selfisolation) and a reduced rate of forming new contacts λ + j = 0.3 • λ + j .Both panels depict results from 10 3 SSATAN-X simulations for each parameter set. A. Error bars represent the mean amount ± standard deviation of the amount of edges eligible for transmission of the virus (i.e.number of contacts between infected and susceptibles and between diagnosed and susceptibles) over time.B. Error bars represent the mean amount ± standard deviation of newly infected individuals over time.Annotations on both panels show the percentage of infected individuals that are diagnosed (D)/(I + D).

Figure 8 .
Figure 8. Run time performance.A. Log-log plot showing the mean (± standard deviation) run time of SSATAN-X (blue error bars) vs. SSA (red error bars) for different population sizes (200, 500 and 1000 individuals) over 10 3 simulations for each parameter set on a single 2x AMD Epyc 7742 64-Core compute node with a base clock of 2.25 Ghz respectively.B. Semi-log plot showing mean (± standard deviation) run time of SSATAN-X (blue error bars) vs. SSA (red error bars) for different parameterizations for a network with 200 individuals.The relation γ γ + (λ − ) 2 , with γ being the transmission rate and λ − being the rate of edge disassembling, can be interpreted as the probability of transmitting an infection before the contact is dissolving.

Figure 9 .
Figure 9.Comparison of simulation results from parallel updating schemes vs. SSATAN-X.A.-C. Central tendency of predictions: Average number of infected, diagnosed and susceptible individuals after 1000 SSATAN-X or parallel updating simulations with different time steps.D.-F.Statistical range of predictions: The dashed lines denote the levels of the average infected, diagnosed and susceptible plus (upper line), or minus (lower line) one standard deviation after 1000 SSATAN-X or parallel updating simulations with different time steps.

Figure 10 .
Figure 10.Kolmogorov-Smirnov test statistics for SSATAN-X and parallel updating in comparison to SSA.The dashed lines show the Kolmogorov-Smirnov Test (KST) statistic of parallel updating (with step sizes dt = 0.0001, 0.1, 0.2 and 0.5), as well as SSATAN-X for the distribution of the number of infected, diagnosed and susceptible individuals.The grey shaded are denotes the baseline KST statistic due to the sampling error (for 10 3 SSA simulations).

simulation of epidemics on adaptive networks bulk contact update leap epidemic change (B L )- bulk contact update leap rejection (no epi. change) acceptance (epi.change)
If the leap is rejected, the sampled values are stored (lines 36-41) and τ is decreased by multiplying it with probability p.The algorithm returns to line 4. Calculate the cumulative sums ς + and ς − for reaction types ξ ∈ {+, −}.Permute M − reactions that eliminate an edge and M + reactions that create an edge and save the order in the vector O 2 for each contact update in O do ).Otherwise, for each reaction type, the number of executions M ξ is sampled (lines 9-19), either from a Poisson distribution (line 12) or a Binomial distribution (line 17).If the leap is accepted (line 21) with condition |E(t + τ )| − |E(t)| ≤ max( |E(t)|, 1) |E (t + τ )| − |E (t)| ≤ max( |E (t)|, 1), (2.4) the time t is increased by τ (line 28), storage variables Q ξ , S ξ and C ξ are updated (lines 22-27) and τ is updated (lines 29-33).The bulk update (line 34) is performed with respect to values of M ξ as outlined in Algorithm 4.1