NoMaSS: A Practical and Scalable Framework for Non-Markovian Gillespie-based Stochastic Simulation

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 the implementattion of a general and scalable framework to simulate non-Markovian stochastic systems with arbitrary inter-event time distribution, mixing delay based methods, MOSAIC and Laplace Gillespie [1].

 

Simulating a non-Markovian system

First, you need to install NoMaSS, or you can use the NoMaSS.py file provided in the repository:

- pip install nomass

Then, you can run a non-Markovian simulation with the toy example below. Other examples, including the three biochemical systems described in the paper: Cell division, differentiation and RNA transcription, are provided in the /NoMaSS/Examples and /REGIR/Biochemical_applications folders.

import nomass

#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 = nomass.Reaction_channel(param,rate=r1, shape_param=alpha1, distribution = 'Gamma', method = 'Delay')
reaction1.reactants = ['A']
reaction1.products = ['B']	
reaction2 = nomass.Reaction_channel(param,rate=r2, shape_param=alpha2, distribution = 'Weibull', method = 'MOSAIC')
reaction2.reactants = ['B']
reaction2.products = ['C','A']
reaction3 = nomass.Reaction_channel(param, rate=r3, shape_param=alpha2, distribution = 'Gamma', method = 'Laplace')
reaction3.reactants = ['C','D']
reaction3.products = []
reaction3 = nomass.Reaction_channel(param, rate=r3)
reaction3.reactants = ['A','B']
reaction3.products = []

#Define the initial population of reactants:
N_init = {'A':300,'B':0,'C':0,'D':0}

#Initialize and run the Gillepsie simulation:
SSA_simul = nomass.Gillespie_simulation(N_init,param)
SSA_simul.reaction_channel_list = [reaction1, reaction2, reaction3]
populations = SSA_simul.run_simulations(param.Tend, reaction_channel_list, verbose = True)

#Plot the results:
SSA_simul.plot_inter_event_time_distribution()
SSA_simul.plot_populations()

The algorithm runs for a few seconds and output the following figures (note that you can disables all printing by passing the argument verbose = False when running the simulation):

The oscillations resulting from the markovian dynamics are clearly visible. If you check carefully, you will notice that the theoretical distributions do not match exactly the simulated distributions, even if you increase the number of simulations. This happens because the entities A and B are reactants of two reaction channels at the same time, and the theoretical distribution only represent the inter-event time distribution that the reaction channel would have if it was the only process interaction with that reactant. In practice, these kind of situations will occur frequently in non-Markovian systems, so do not worry if the simulated and theoretical distributions do not match exactly. The accuracy of each method was rigourously demonstrated in [1] and [2].

State dependant waiting time distributions

You can continuously update IED parameters during the simulation via the custom_rates_update callback (called at each iteration / after each event, depending on your setup). This is useful for state-dependent kinetics, where reaction parameters depend on the current system state (e.g., population sizes, resource levels, feedback loops).

Important: delay-based methods (DelaySSA) can become inaccurate when parameters change over time, because scheduled delays are drawn at initiation and cannot be revised once they are placed in the queue. This effectively “freezes” the distribution parameters until the event fires, which is not correct for rapidly changing systems. For time-varying or state-dependent distributions, prefer Laplace Gillespie or MOSAIC, which remain exact under parameter updates (they resample / evaluate rates based on the current state).

import nomass
[...]
def custom_rates_update(SSA_simul):
    #A custom function to update the IED parameters of each channels   
    
    NA = len(SSA_simul.reactant_list['A'])
	NB = len(SSA_simul.reactant_list['B'])
    N = NA + NB
    
    SSA_simul.reaction_channel_list[2].rate = r1 * (1 - N/100)
    #SSA_simul.reaction_channel_list[2].shape_param = 0.5 # You can also adjust the shape param as you whish, although most models will typically only involve changing the rate.
	
[...]
populations = SSA_simul.run_simulations(param.Tend, reaction_channel_list, update_function = custom_rates_update, verbose = True)

Implemented distributions

NoMaSS implement three main class of methods: MOSAIC [2], Laplace Gillespie [3] and DelaySSA [4]

PDF name	MOSAIC [2]	Laplace [3]	DelaySSA [4]
Exponential	✅	✅	✅
Gamma	✅ (α ≥ 1)	✅ (α ≤ 1)	✅
Weibull	✅ (α ≥ 1)	✅ (α ≤ 1)	✅
Gaussian / Normal	✅	❌	✅
Lognormal	✅	❌	✅
Cauchy	✅	❌	✅
Power-law (Pareto)	✅	❌	✅
Power-law (Lomax)	✅	✅	✅

With the current implementation, each distribution is parameterized by a nominal rate λ (inverse mean waiting time — or median when the mean is undefined) and a shape parameter α, as detailed below [Full IED parametrizations are provided in reference [1] Table S1.]:

  Exponential (https://en.wikipedia.org/wiki/Exponential_distribution):
      - rate: 1/mean
      - shape: None
  
  Normal / Gaussian (https://en.wikipedia.org/wiki/Normal_distribution):
      - rate: 1/mean
      - shape: std/mean
  
  LogNormal (https://en.wikipedia.org/wiki/Log-normal_distribution):
      - rate: 1/mean
      - shape: std/mean
      
  Gamma (https://en.wikipedia.org/wiki/Gamma_distribution):
      - rate: 1/mean
      - shape: α
      
  Weibull (https://en.wikipedia.org/wiki/Weibull_distribution):
      - rate: 1/mean
      - shape: k
      
  Cauchy (https://en.wikipedia.org/wiki/Cauchy_distribution):
      - rate: 1/median
      - shape: std/median

  Pareto type 1 (https://en.wikipedia.org/wiki/Pareto_distribution):
      - rate: 1/median
      - shape: α

  Lomax / Pareto type 2 (https://en.wikipedia.org/wiki/Lomax_distribution):
      - rate: 1/median
      - shape: α

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.

Feel free to drop me an email if you have any question, I will be happy to help !

Customizing NoMaSS for your system

The NoMaSS framework offer countless possibilities and highly customizable models. However, with the current implementation, reactions propensities are always proportional to the number of reactant. For example, the reaction (A+B -> C) will have a propensity of a = A x B x r. In some models, you might want to implement more complex formula for the reaction propensities, (such as for example a = A x B x r / D, where D is a parameter that evolves with the system). To do so, you can directly modify the NoMaSS/compute_propensities function according to your need. Likewise, you might want to set up specific rejection rules if your reactants have some individual properties, and modify them appropriatly. To do so, first define your reactant properties in the NoMaSS/Reactant class, and then define your reaction specific rules in NoMaSS/perform_reaction.

Feel free to drop me an email if you are not sure how to do it, I will be happy to help !

References

[1] Pélissier, A, et al. "Gillespie-Based Methods for Non-Markovian Stochastic Simulation: A Practical Guide". Under Review (2026)

[2] Pélissier, A, et al. "Unifying non-Markovian Dynamics and Agent Heterogeneity in Scalable Stochastic Networks". Nature Communication (2026)

[3] Masuda, Naoki, and Luis EC Rocha. "A Gillespie algorithm for non-Markovian stochastic processes." Siam Review 60.1 (2018): 95-115.

[4] David F Anderson. “A modified next reaction method for simulating chemical systems with time dependent propensities and delays”. The Journal of chemical physics 127.21 (2007).