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 and efficient algorithm like the Gillespie algorithm for Markovian processes, with the difference being 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 nonexponential interevent 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 because it can be implemented very easily in cases (such as for delays coupled to the evolution of the system) in which other algorithms do not work or need adapted versions that are less efficient in computational terms.

Figure 1

Testing the algorithm. Survival probability for three processes generated by the nMGA with parameters $(\alpha =1,\mu =1)$, $(\alpha =1.5,\mu =1)$, and $(\alpha =0.5,\mu =0.1)$. Solid lines are given by Eq. (11) with the same parameters. In (a) the total number of processes is $N={10}^3$, whereas in (b) it is only $N=20$.

Figure 2

Non-Markovian epidemic spreading. Top: Two possible ways to generate a new active link in the SIS model. Rules 1 and 2 specify the elapsed time for the newborn active link. Rule 1 sets this time always to zero, whereas rule 2 assigns to the newborn active link the elapsed time of individual A. Bottom: Prevalence of the epidemics at the steady state as a function of the effective infection rate in Erdös-Rényi networks of size $N={10}^4$ and average degree $\langle k\rangle =5$. Solid lines stand for non-Markovian dynamics with rule 1 implemented, whereas dashed lines correspond to rule 2. The red solid curve corresponds to the Markovian case. In all cases, recovery events are exponentially distributed.

Figure 3

Cooperative non-Markovian epidemic spreading. Prevalence of the epidemics at the steady state as a function of the effective infection rate for the same network as in Fig. 2 in the case of cooperative infections given by Eq. (13) for different values of $\alpha$ and $\sigma$. Rule 1 is used in all cases.

Figure 4

Performance of the nMGA and modified Anderson's algorithms for the gene regulation model Eq. (14) with time delays coupled to the evolution of the system. Top: Scaling of computational times required to simulate ${10}^4$ physical time units, starting with $n_{M*}=0$ and the steady-state values of $n_p$ and $n_M$ for $\alpha =5$. Bottom: Autocorrelation function, normalized by system size, for the temporal evolution of the number of proteins in the stationary state as a function of the shape parameter $\alpha$ of the time-delay distribution, for $N=100$. In both cases, the delays follow a Weibull distribution with the parameters detailed in the main text.