Source-linked AI summary
Simulating non-Markovian stochastic processes
Marian Boguna, Luis F. Lafuerza, Raul Toral, M. Angeles Serrano
TL;DR
Non-exponential inter-event times introduce memory effects that make many stochastic processes non-Markovian and difficult to analyze. The paper provides an exact framework and a Gillespie-like algorithm for simulating such processes, then applies it to epidemic spreading and delayed biochemical reactions. Its results show that subtle differences in non-Markovian rules can substantially affect complex-system behavior, and that the algorithm applies to delay-coupled cases where existing methods cannot be straightforwardly used.
Problem
Non-exponential inter-event times occur in many settings, creating memory effects that make processes non-Markovian and difficult to analyze.
Source: §I · p. 1Method
The paper gives an exact framework and a large-system, Gillespie-like algorithm for simulating discrete stochastic processes with arbitrary inter-event distributions.
Source: §I · p. 1, §VI · p. 7Results
Subtle differences in non-Markovian dynamical rules have drastic effects on complex systems’ global behavior, and the nMGA applies to delay-coupled biochemical reactions where prior methods cannot be straightforwardly applied.
Source: §V · p. 6, §VI · p. 7Takeaways & Limitations
Correctly modeling non-Markovian events matters because minimal variations in their dynamics can have drastic effects on global behavior.
Source: §VI · p. 7
Abstract
from arXiv · showhide
We present a simple and general framework to simulate statistically correct realizations of a system of non-Markovian discrete stochastic processes. We give the exact analytical solution and a practical an efficient algorithm alike the Gillespie algorithm for Markovian processes, with the difference that now the occurrence rates of the events depend on the time elapsed since the event last took place. We use our non-Markovian generalized Gillespie stochastic simulation methodology to investigate the effects of non-exponential inter-event time distributions in the susceptible-infected-susceptible model of epidemic spreading. Strikingly, our results unveil the drastic effects that very subtle differences in the modeling of non-Markovian processes have on the global behavior of complex systems, with important implications for their understanding and prediction. We also assess our generalized Gillespie algorithm on a system of biochemical reactions with time delays. As compared to other existing methods, we find that the generalized Gillespie algorithm is the most general as it can be implemented very easily in cases, like for delays coupled to the evolution of the system, where other algorithms do not work or need adapted versions, less efficient in computational terms.
I. INTRODUCTION
Non-exponential inter-event times create memory effects that make many discrete stochastic processes non-Markovian, and such processes are difficult to analyze. The paper develops an exact framework and an efficient Gillespie-like algorithm for simulating them, then applies it to epidemic spreading and delayed biochemical reactions.
I. INTRODUCTION
- Non-exponential inter-event times produce memory effects and non-Markovian dynamics, which are difficult to understand analytically and often studied through numerical simulation.Reported examples span reliability, queuing, human activity, temporal networks, and biochemical reactions with delays.Source: §I · p. 1
- The paper develops a framework for statistically correct simulations of discrete processes with arbitrary inter-event time distributions.Processes may stochastically create or annihilate other processes and depend on the system’s current state.Source: §I · p. 1
- An exact solution and an efficient, simple simulation algorithm extend the practical spirit of Gillespie’s method to these processes.The paper applies the algorithm to susceptible-infected-susceptible epidemic spreading and coupled biochemical reactions with time delays.Source: §I · p. 1
- Subtle differences in non-Markovian dynamical rules can have important effects on the global behavior of complex systems.Source: §I · p. 1
II. MARKOVIAN STOCHASTIC SIMULATION: THE GILLESPIE ALGORITHM
The Gillespie algorithm simulates coupled discrete Markovian processes by repeatedly drawing the next event time and selecting which process occurs. Changes to the system can also alter the rates and active processes used in subsequent steps.
II. MARKOVIAN STOCHASTIC SIMULATION: THE GILLESPIE ALGORITHM
- At each step, the algorithm draws the time until the next event and selects the process that occurs.For independent processes with rates λ_i, the method uses a waiting-time distribution and a process-selection probability.Source: §II · p. 1
- An event can change the system state, the rates, or which processes are active for the next iteration.For example, a reaction changes molecule counts and therefore the rates of reactions involving those species.Source: §II · p. 2
- The method can simulate reaction-like systems where the realization stochastically changes the number of processes, including dynamics with absorbing states.Source: §II · p. 2
III. NON-MARKOVIAN STOCHASTIC SIMULATION
The non-Markovian simulation conditions each process’s next-event probability on the time elapsed since its last event. It samples the next event time and process, updates elapsed times and system state, and initializes elapsed times according to the application.
III. NON-MARKOVIAN STOCHASTIC SIMULATION
- For each process, the probability density of its next event depends on its inter-event distribution and the time elapsed since its previous event.The section distinguishes elapsed time t_i from future waiting time τ.Source: §III · p. 2
- The joint probability for the next event accounts for the chosen process occurring while every other process does not occur during the same interval.Conditional survival probabilities are formed from the survival functions and elapsed times of the processes.Source: §III · p. 2
- The algorithm draws a waiting time from the cumulative distribution, selects a process from the resulting discrete distribution, and updates elapsed times.Elapsed times advance by τ for processes that did not fire, while the selected process’s elapsed time resets to zero.Source: §III · p. 2
- After an event, the algorithm updates the system state and active processes, assigning elapsed times to any newly activated processes.Source: §III · p. 2
- Elapsed times can start at zero or be initialized from the steady-state age distribution Ψ_i(t)/⟨τ_i⟩.The latter option assumes the system is already in steady state.Source: §III · p. 2
A. Generalized Gillespie algorithm
The large-system approximation yields a Gillespie-like algorithm using elapsed-time-dependent rates. Tests against the theoretical survival distribution show accurate inter-event times, including for a small number of processes, with a small deviation when α_i > 1.
A. Generalized Gillespie algorithm
- When N is large, the waiting-time distribution can be simplified by expanding the survival probability near τ = 0.The approximation becomes exact in the limit N → ∞.Source: A. Generalized Gillespie algorithm · pp. 2–3
- The approximation omits the process that fired most recently, making consecutive selections of the same process impossible.The authors state that this event is negligible for large N, while the exclusion avoids potential divergence when limτ→0+ ψ_i(τ) = ∞.Source: A. Generalized Gillespie algorithm · p. 3
- The resulting nMGA recovers the classical Gillespie algorithm for constant rates and otherwise uses rates that depend on elapsed times.Source: A. Generalized Gillespie algorithm · p. 3
- The test mixes Weibull renewal processes with scale and shape parameters drawn uniformly from µ_i ∈ (0.1, 1) and α_i ∈ (0.5, 1.5).The simulations compare inter-event survival probabilities for selected processes against their corresponding theoretical distributions.Source: A. Generalized Gillespie algorithm · p. 3
- For N = 20, the nMGA reproduces inter-event times accurately, with a small deviation for processes with α_i > 1.The figure compares simulated survival probabilities with the theoretical distribution; its panels use N = 10^3 and N = 20.Source: A. Generalized Gillespie algorithm · p. 3, Figure 1 · p. 3
IV. EPIDEMIC SPREADING. THE SIS MODEL AS A CASE STUDY
The SIS model represents individuals as susceptible or infected: infected individuals recover, while susceptible individuals can be infected by infected neighbors. On Erdős–Rényi networks, the model has an absorbing-to-endemic phase transition at a topology-dependent effective infection rate.
IV. EPIDEMIC SPREADING. THE SIS MODEL AS A CASE STUDY
- In the SIS model, infected individuals recover spontaneously, while susceptible individuals are infected at a rate proportional to their infected neighbors.Source: §IV · p. 4
- The dynamics transitions from an absorbing healthy phase, where outbreaks disappear exponentially fast, to an endemic phase with sustained activity.Source: §IV · p. 4
- On Erdős–Rényi networks, the critical effective infection rate depends on network topology and is approximately 1/⟨k⟩ for large networks.Source: §IV · p. 4
A. Independent infections
With independent active-link infection processes, non-Markovian timing and the rule for initializing new links both affect epidemic prevalence. The effect of non-Poisson recovery is smaller than that of non-Markovian infection timing in the cases considered.
A. Independent infections
- Each active link is an independent infection process, and a susceptible with multiple infected neighbors is infected when the first link's event fires.In non-Markovian dynamics, infection depends on each link's elapsed time, not just the number of infected neighbors.Source: A. Independent infections · p. 4
- Non-Markovian infection timing shifts the epidemic critical point, and for the same Weibull parameter α, prevalence under rule 2 is above rule 1 when α > 1 and below it when α < 1.Rule 1 resets a new link's elapsed time to zero; rule 2 assigns it the infected individual's elapsed infection time.Source: A. Independent infections · pp. 4–5, Figure 2 · p. 4
- For Weibull infection times, the mean time until the next event increases with elapsed time when α < 1 and decreases when α > 1, explaining rule 2's longer or shorter infection times relative to rule 1.Source: A. Independent infections · p. 5
- Non-Poisson recovery has a much smaller effect than non-Markovian infection: rule 1 shows no noticeable difference from the Markovian case, while rule 2 differs only for very heterogeneous recovery times.Source: A. Independent infections · p. 5
B. Cooperative infections
The cooperative-infection extension allows a susceptible individual's total infection rate to depend nonlinearly on the active links. Nonlinearity changes prevalence but does not shift the critical point in the cases described.
B. Cooperative infections
- The independent-link assumption is reasonable under rule 2, but rule 1 leaves room for active links to cooperate nonlinearly.Source: B. Cooperative infections · p. 5
- The total infection rate combines active-link contributions through a nonlinear exponent σ; σ = 1 recovers independent infections.The rate uses the adjacency matrix, infected-neighbor states, and the duration of each active link.Source: B. Cooperative infections · p. 5
- Nonlinear infection increases prevalence when σ > 1 and decreases it when σ < 1, without changing the critical point.Near the critical point, few infected neighbors make the dynamics effectively independent.Source: B. Cooperative infections · p. 5
V. ASSESSMENT OF THE GENERALIZED GILLESPIE ALGORITHM ON A SYSTEM OF BIOCHEMICAL REACTIONS WITH DELAYS
The section tests the generalized Gillespie algorithm (nMGA) on delayed gene regulation, including delays whose distribution changes with system state. In this setting, nMGA is faster than the adapted Anderson method, while the model’s protein autocorrelation varies with the delay-distribution shape.
V. ASSESSMENT OF THE GENERALIZED GILLESPIE ALGORITHM ON A SYSTEM OF BIOCHEMICAL REACTIONS WITH DELAYS
- When delays are independent of system evolution, nMGA does not necessarily outperform prior methods; when they are state-coupled, earlier methods cannot be applied straightforwardly and adapted Anderson is slower than nMGA.A separate numerical-integration approach may suit very large systems, but resolving the highest frequencies raises its computing time, and it does not directly apply to networked populations.Source: §V · pp. 6–7
- The model represents delayed auto-inhibition: protein-dependent mRNA transcription initiation is followed by delayed completion, while mRNA and protein degradation and protein translation are Poissonian.The section also motivates state-coupled delays by finite transcription resources, with delay time affected by mRNA being transcribed or already present.Source: §V · p. 6
- Both algorithms’ computational times scale as N^2, while nMGA is ≈3.7 times faster for α = 0.5 and ≈9.5 times faster for α = 5.These comparisons concern delays coupled to system evolution; the adapted Anderson method performs relatively well because the Weibull cumulative distribution has an elementary-function expression.Source: Figure 4 · p. 6, §V · p. 6
- As α rises from 0.5 to 100, the increasingly peaked delay distribution accompanies more pronounced decaying oscillations in protein autocorrelation, which eventually disappear as delays become burstier.The autocorrelation is for protein number in the stationary state and is normalized by system size; the wave cycle is marked by the average delay.Source: §V · pp. 6–7, Figure 4 · p. 6
VI. CONCLUSIONS
The paper introduces an exact framework for simulating non-Markovian discrete stochastic processes and a simpler large-system algorithm with Gillespie-like operation. Its results show that small changes in non-Markovian event modeling can have drastic effects on system-level behavior.
VI. CONCLUSIONS
- The framework generates statistically correct realizations of non-Markovian discrete stochastic processes, with a simple simulation algorithm as its large-system approximation.The algorithm resembles Gillespie’s method, except event rates depend on elapsed time since the event last occurred.Source: §VI · p. 7
- The algorithm is not always the fastest, but it is especially general for delays coupled to system evolution, where other methods may not apply directly or need adaptation.Source: §VI · p. 7
- Subtle differences in non-Markovian modeling can have drastic effects on global system behavior, making correct modeling important for understanding and prediction.Source: §VI · p. 7
Appendix A: Simulation details for the SIS dynamics
The appendices describe event-list updates for SIS simulations and event selection and elapsed-time tracking for delayed biochemical reactions.
Appendix A: Simulation details for the SIS dynamics
- For SIS dynamics, the simulation selects recovery or infection from lists of infected nodes and active links, draws an exponential time increment, and updates the lists after each event.New infections enter with elapsed time zero; recovery and infection alter which links are active, and remaining elapsed times can be computed from birth times.Source: Appendix A · pp. 7–8
- For SIS cooperative infections, the procedure adds a list of susceptible nodes with rates from Eq. (13), then chooses events from infected-node and susceptible-node lists.Source: Appendix A · p. 8
Appendix B: Simulations details for delayed reactions
- For delayed reactions, the simulation tracks protein and mRNA counts plus in-progress transcription events with elapsed times and state-dependent completion rates, then selects among five event types.Newly initiated transcription is excluded from event selection at elapsed time zero; after a time increment, ongoing transcription elapsed times and affected molecule counts and lists are updated.Source: Appendix B · p. 8
Appendix C: Adaptation of Anderson’s modified next reaction method for general non-exponential inter-event times
The appendix adapts Anderson’s modified next reaction method to event rates that depend on elapsed time since the reaction last fired. It uses internal times to determine candidate firing times and iterates after each event.
Appendix C: Adaptation of Anderson’s modified next reaction method for general non-exponential inter-event times
- For reaction i, the state changes by its stoichiometric increment when it fires, at a rate depending on current state and elapsed time since its last firing.Source: Appendix C · p. 8
- Independent unit-rate Poisson processes and exponential internal times encode candidate reaction firings, with firing times obtained by solving the integrated rate relation.The internal times are independent of system state; after setting them, firing times and state changes follow from the relation.Source: Appendix C · p. 8
- To find the next event, each reaction’s candidate time is solved using the current state, and the earliest candidate fires; elapsed internal time then supports the next iteration.Source: Appendix C · pp. 8–9