An algorithm for non-Markovian stochastic simulation
Project description
REGIR: A Scalable Gillespie Algorithm for non-Markovian Stochastic Simulations
Refer to the github folder for a full description (https://github.com/Aurelien-Pelissier/REGIR). Discrete stochastic processes are widespread in both nature and human-made systems, with applications across physics, biochemistry, epidemiology, social patterns and finance, just to name a few. In the majority of these systems, the dynamics cannot be properly described with memoryless (or Markovian) interactions, and thus require the use of numerical tools for analyzing these non-Markovian dynamics. This repository contains an implementattion of a general and scalable framework to simulate non-Markovian stochastic systems with arbitrary inter-event time distribution and accuracy. The algorithm is referred to as the Rejection Gillespie algorithm for non-Markovian Reactions (REGIR) [1].
Simulating a non-Markovian system
First, you need to install REGIR, or you can use the REGIR.py
file provided in the repository:
- pip install REGIR
Then, you can directly run a non-Markovian simulation with this toy example (other examples, including the three biochemical systems described in the paper: Cell division, differentiation and RNA transcription, are provided in the /REGIR/Examples
folder.):
import REGIR as gil
#Set the simulation parameters:
class param:
Tend = 10 #Length of the simulation
N_simulations = 20 #The simulation results should be averaged over many trials
unit = 'h' #Unit of time (used for plotting only)
timepoints = 100 #Number of timepoints to record (used for plotting only)
r1 = 1
r2 = 4
r3 = 0.03
alpha1 = 20
alpha2 = 5
#Define the reaction chanels:
reaction1 = gil.Reaction_channel(param,rate=r1, shape_param=alpha1, distribution = 'Gamma')
reaction1.reactants = ['A']
reaction1.products = ['B']
reaction2 = gil.Reaction_channel(param,rate=r2, shape_param=alpha2, distribution = 'Weibull')
reaction2.reactants = ['B']
reaction2.products = ['C','A']
reaction3 = gil.Reaction_channel(param,rate=r3)
reaction3.reactants = ['A','B']
reaction3.products = []
#Define the initial population of reactants:
N_init = dict()
N_init['A'] = 300
N_init['B'] = 0
N_init['C'] = 0
#Initialize and run the Gillepsie simulation:
reaction_channel_list = [reaction1, reaction2, reaction3]
G_simul = gil.Gillespie_simulation(N_init,param)
G_simul.reaction_channel_list = reaction_channel_list
G_simul.run_simulations(param.Tend, verbose = True)
G_simul.plot_inter_event_time_distribution()
G_simul.plot_populations()
The algorithm runs for a few seconds and output the following figures (note that you can disables all printing and plotting by passing the argument verbose = False
when running the simulation)
Implemented distributions
With the current implementation, each available distribution are characterised by their rate and a shape parameter as follow:
Exponential:
- rate: 1/mean
- shape parameter: None
Normal:
- rate: 1/mean
- shape: std/mean
LogNormal:
- rate: 1/mean
- shape: std/mean
Gamma:
- rate: 1/mean
- shape: α >= 1 (https://en.wikipedia.org/wiki/Gamma_distribution)
Weibull:
- rate: 1/mean
- shape: k >= 1 (https://en.wikipedia.org/wiki/Weibull_distribution)
Cauchy:
- rate: 1/median
- shape: γ (https://en.wikipedia.org/wiki/Cauchy_distribution)
Keep in mind that non-Markovian simulations are only available for reaction channels with a single reactant, as the definition of inter-event time distribution is ambigious for channels with multiple reactants. If a channel is defined without or with more than one reactant, it will be considered as a Poisson process. Also, note that monotolically decreasing distributions, such as Weibull (k < 1), gamma (α < 1) or power laws, are not available in the current implementation of this repository, as these can be more elegantly and efficiently simulated with the Laplace Gillespie algorithm [2].
Feel free to drop me an email if you have interest in me adding the Laplace Gillespie or any other relevant distributions to this implementation.
References
[1] Pelissier, A, Phan, M, et al. "Practical and scalable simulations of non-Markovian stochastic processes". Proceedings of the National Academy of Sciences (2022)
[2] Masuda, Naoki, and Luis EC Rocha. "A Gillespie algorithm for non-Markovian stochastic processes." Siam Review 60.1 (2018): 95-115.
Project details
Release history Release notifications | RSS feed
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.