UC-Davis-molecular-computing / ppsim

Python package for simulating population protocols
MIT License
10 stars 5 forks source link

allow timescale to be configurable to match between different methods of specifying transition rules #38

Closed dave-doty closed 2 years ago

dave-doty commented 2 years ago

The following code runs the same protocol, with one transition (A,A) → (B,B), using the CRN way of specifying transitions and the dictionary way (the former always uses continuous time until we fix issue #37, so we set the latter to also use continuous time):

import ppsim as pp

a,b = pp.species('A B')
cfg = {a: 10**6}
rxns = [2*a >> 2*b]
transitions = { (a,a): (b,b) }
sim_crn = pp.Simulation(cfg, rxns)
sim_pp = pp.Simulation(cfg, transitions, continuous_time=True)
sim_crn.run(1, 0.01)
sim_pp.run(1, 0.01)
sim_crn.history.plot()
sim_pp.history.plot()

However, the timescale is about factor-2 different between them:

image

See https://colab.research.google.com/drive/1jXSnNp8mgReIq74hYkZ1UW6o2jtwEqzm#scrollTo=Zi522V7PnbPS to run the code.

The issue appears to be that for the timescales to match, one would have to define 1 unit of time in a population protocol to be (n-1)/2 interactions, but then it would not match the standard literature definition of n interactions = 1 unit of time. Alternately, one could apply a factor ≈2 speedup to the CRN model, but then it would not match standard Gillespie kinetics.

A good compromise might be to either do the former (scale time for non-CRN protocols to be 1 unit of time = (n-1)/2 interactions) by default, or just leave it as is, but give a parameter to the Simulation constructor that lets the user re-define the time scaling if they want. Not sure how fancy to get with that parameter. It could be just a Boolean indicating whether to use (n-1)/2 or n, or it could be something like a function that takes as input a number of interactions i and the population size n and outputs a float indicating how much time it represents.

dave-doty commented 2 years ago

Closing, since I forgot this is already configurable by the user by setting the parameter steps_per_time_unit in the Simulation constructor.