Simulating non-Markovian stochastic processes

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.Comment: Improvement of the algorithm, new results, and a major reorganization of the paper thanks to our coauthors L. Lafuerza and R. Tora

Temporal Gillespie algorithm: Fast simulation of contagion processes on time-varying networks

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

A gillespie algorithm for non-markovian stochastic processes

The Gillespie algorithm provides statistically exact methods for simulating stochastic dynamics modeled as interacting sequences of discrete events including systems of biochemical reactions or earthquake occurrences, networks of queuing processes or spiking neurons, and epidemic and opinion formation processes on social networks. Empirically, the inter-event times of various phenomena obey long-tailed distributions. The Gillespie algorithm and its variants either assume Poisson processes (i.e., exponentially distributed inter-event times), use particular functions for time courses of the event rate, or work for non-Poissonian renewal processes, including the case of long-tailed distributions of inter-event times, but at a high computational cost. In the present study, we propose an innovative Gillespie algorithm for renewal processes on the basis of the Laplace transform. The algorithm makes use of the fact that a class of point processes is represented as a mixture of Poisson processes with different event rates. The method is applicable to multivariate renewal processes whose survival function of inter-event times is completely monotone. It is an exact algorithm and works faster than a recently proposed Gillespie algorithm for general renewal processes, which is exact only in the limit of infinitely many processes. We also propose a method to generate sequences of event times with a tunable amount of positive correlation between inter-event times. We demonstrate our algorithm with exact simulations of epidemic processes on networks, finding that a realistic amount of positive correlation in inter-event times only slightly affects the epidemic dynamics

Practical and scalable simulations of non-Markovian stochastic processes

Discrete stochastic processes are widespread in natural systems with many applications across physics, biochemistry, epidemiology, sociology, and finance. While analytic solutions often cannot be derived, existing simulation frameworks can generate stochastic trajectories compatible with the dynamical laws underlying the random phenomena. However, most simulation algorithms assume the system dynamics are memoryless (Markovian assumption), under which assumption, future occurrences only depend on the present state of the system. Mathematically, the Markovian assumption models inter-event times as exponentially distributed variables, which enables the exact simulation of stochastic trajectories using the seminal Gillespie algorithm. Unfortunately, the majority of stochastic systems exhibit properties of memory, an inherently non-Markovian attribute. Non-Markovian systems are notoriously difficult to investigate analytically, and existing numerical methods are computationally costly or only applicable under strong simplifying assumptions, often not compatible with empirical observations.To address these challenges, we have developed the Rejection-based Gillespie algorithm for non-Markovian Reactions (REGIR), a general and scalable framework to simulate non-Markovian stochastic systems with arbitrary inter-event time distributions. REGIR can achieve arbitrary user-defined accuracy while maintaining the same asymptotic computational complexity as the Gillespie algorithm. We illustrate REGIR's modeling capabilities in three important biochemical systems, namely microbial growth dynamics, stem cell differentiation, and RNA transcription. In all three cases, REGIR efficiently models the underlying stochastic processes and demonstrates its utility to accurately investigate complex non-Markovian systems. The algorithm is implemented as a python library REGIR

Stochastic Wilson-Cowan models of neuronal network dynamics with memory and delay

We consider a simple Markovian class of the stochastic Wilson-Cowan type models of neuronal network dynamics, which incorporates stochastic delay caused by the existence of a refractory period of neurons. From the point of view of the dynamics of the individual elements, we are dealing with a network of non-Markovian stochastic two-state oscillators with memory which are coupled globally in a mean-field fashion. This interrelation of a higher-dimensional Markovian and lower-dimensional non-Markovian dynamics is discussed in its relevance to the general problem of the network dynamics of complex elements possessing memory. The simplest model of this class is provided by a three-state Markovian neuron with one refractory state, which causes firing delay with an exponentially decaying memory within the two-state reduced model. This basic model is used to study critical avalanche dynamics (the noise sustained criticality) in a balanced feedforward network consisting of the excitatory and inhibitory neurons. Such avalanches emerge due to the network size dependent noise (mesoscopic noise). Numerical simulations reveal an intermediate power law in the distribution of avalanche sizes with the critical exponent around -1.16. We show that this power law is robust upon a variation of the refractory time over several orders of magnitude. However, the avalanche time distribution is biexponential. It does not reflect any genuine power law dependence

Stochastic processes with distributed delays: chemical Langevin equation and linear-noise approximation

We develop a systematic approach to the linear-noise approximation for stochastic reaction systems with distributed delays. Unlike most existing work our formalism does not rely on a master equation, instead it is based upon a dynamical generating functional describing the probability measure over all possible paths of the dynamics. We derive general expressions for the chemical Langevin equation for a broad class of non-Markovian systems with distributed delay. Exemplars of a model of gene regulation with delayed auto-inhibition and a model of epidemic spread with delayed recovery provide evidence of the applicability of our results.Comment: 21 pages, 7 figure

Optimized Gillespie algorithms for the simulation of Markovian epidemic processes on large and heterogeneous networks

Numerical simulation of continuous-time Markovian processes is an essential and widely applied tool in the investigation of epidemic spreading on complex networks. Due to the high heterogeneity of the connectivity structure through which epidemics is transmitted, efficient and accurate implementations of generic epidemic processes are not trivial and deviations from statistically exact prescriptions can lead to uncontrolled biases. Based on the Gillespie algorithm (GA), in which only steps that change the state are considered, we develop numerical recipes and describe their computer implementations for statistically exact and computationally efficient simulations of generic Markovian epidemic processes aiming at highly heterogeneous and large networks. The central point of the recipes investigated here is to include phantom processes, that do not change the states but do count for time increments. We compare the efficiencies for the susceptible-infected-susceptible, contact process and susceptible-infected-recovered models, that are particular cases of a generic model considered here. We numerically confirm that the simulation outcomes of the optimized algorithms are statistically indistinguishable from the original GA and can be several orders of magnitude more efficient.Comment: 12 pages, 9 figure

Markovian Dynamics on Complex Reaction Networks

Complex networks, comprised of individual elements that interact with each other through reaction channels, are ubiquitous across many scientific and engineering disciplines. Examples include biochemical, pharmacokinetic, epidemiological, ecological, social, neural, and multi-agent networks. A common approach to modeling such networks is by a master equation that governs the dynamic evolution of the joint probability mass function of the underling population process and naturally leads to Markovian dynamics for such process. Due however to the nonlinear nature of most reactions, the computation and analysis of the resulting stochastic population dynamics is a difficult task. This review article provides a coherent and comprehensive coverage of recently developed approaches and methods to tackle this problem. After reviewing a general framework for modeling Markovian reaction networks and giving specific examples, the authors present numerical and computational techniques capable of evaluating or approximating the solution of the master equation, discuss a recently developed approach for studying the stationary behavior of Markovian reaction networks using a potential energy landscape perspective, and provide an introduction to the emerging theory of thermodynamic analysis of such networks. Three representative problems of opinion formation, transcription regulation, and neural network dynamics are used as illustrative examples.Comment: 52 pages, 11 figures, for freely available MATLAB software, see http://www.cis.jhu.edu/~goutsias/CSS%20lab/software.htm