Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Unexpected near-duplications in stochastic simulation output with different initial conditions. #117

Open
sclamons opened this issue Aug 14, 2021 · 0 comments

Comments

@sclamons
Copy link
Collaborator

I have a simulation of a very simple birth-death process (DNA replicating uncontrolledly) where I'm running the same model through bioscrape with different initial conditions. When I run these, I usually (for many random seeds) see that at least two traces are, after some short time, identical but shifted in time. It looks like somehow, it's very easy to get into a state where both the RNG and the chemical system are in identical states... which seems like it should, instead, be hard.

Importantly, I set the same random seed at the start of each simulation. If I set the same seed at the start of the first simulation and no others, then I get totally uncorrelated-looking traces, as expected.

Code to replicate (in a Jupyter notebook/lab):

import matplotlib.pyplot as plt
import bioscrape as bs
import numpy as np
%matplotlib inline

rxns = [
    (['DNA'], ['DNA', 'DNA'], 'massaction', {'k': 1}),
    (['DNA'], [], 'massaction', {'k': 1})
]

inits = np.linspace(1, 100, 10)
ts = np.linspace(0, 300, 1000)
stoch_results = []
for i, init_DNA in enumerate(inits):
    seed = 42339
    np.random.seed(seed)
    bs.random.py_seed_random(seed)
    trivial_model = bs.types.Model(species = ["DNA"], 
                               parameters = {}, 
                               reactions = rxns, 
                               initial_condition_dict = {"DNA": init_DNA})
    stoch_results.append(bs.simulator.py_simulate_model(ts, Model = trivial_model, stochastic = True))

for i in range(len(inits)):
    plt.plot(ts, stoch_results[i]["DNA"], 
              lw = 1, ls = "-", label = f"#{i}")
plt.title("DNA replicating with trivial mechanism", fontsize = 16)
plt.xlabel("Time (nondimensionalized)", fontsize = 16)
plt.ylabel("# DNA", fontsize = 16)
plt.legend(loc='upper right')
plt.tight_layout()
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant