Is there any way to get "exact" times of reaction events?

[jonstrutz11]

Maybe this isn't possible due to the implementation details, but is it possible to get a results set where each time is a specific reaction event?

So instead of something like (current output): times:      0, 1, 2, 3, 4, 5 reactant: 10, 9, 7, 4, 2, 0

You could get something like: times:      0, 0.8312, 1.523, 1.852, 2.564, 2.588, 2.911, 3.581, 3.661, 4.111, 4.232 reactant: 10,         9,       8,        7,        6,       5,        4,        3,        2,       1,       0

Or is there a way to "track" whenever a specific species changes by 1?

So far, I tried setting the increment to a very low value, but this is (probably) computationally quite wasteful since during most of these very small timesteps nothing changes. But it does give me relatively accurate values for what I want. Just wondering if there is a faster or more accurate way to do this.

[seanebum]

Hi @jonstrutz11 ,

This is not a something we support, and sampling at small increments as you mentioned is the closest way to achieve this level of granularity. Do you mind if I ask for more information regarding what you are trying to do?

[jonstrutz11]

Hi @seanebum. I am trying to use a KMC model to model a polymerization reaction where there are multiple monomers. (e.g. monomer A, B, and C). So I'm trying to build a generative model that uses kinetic parameters from expt data/parameter inference (in a continuum version of the model) and simulates the generation of polymer chains (formed of A, B, and C) over time. My plan to do this is to run a KMC simulation, get the times of every incorporation event, and then for each event randomly choose a polymer chain to incorporate that monomer. Ideally, I'd be able to tailor the kinetics to not only the monomer identity but also the identity of the last few monomers in the polymer chain it adds to, as those effect the kinetcs as well.

I am also interested in capturing the distribution of times between monomer incorporation events and how that changes based on things like monomer identity, the identity of the previous monomer added to the reactant polymer chain, etc..

Hopefully that isn't too confusing. But I think knowing when every specific monomer incorporation event happens is what I need for this.

[seanebum]

Hi @jonstrutz11, Sorry for the delay! After discussing your model and needs with the team, I do not think GillesPy2 will work directly for what you are attempting. That being said, we would love to collaborate with you on this project if you are interested!

[Princeofthebow]

Hi - I have such a request as well. Somewhere under the hood (for a simple SSA Gillepsie run) the time stamps should be computed? Saving them helps in computing the distribution by weighting the events per time of occurence.

[briandrawert]

For the Numpy solver of GillesPy2, the event times are calculated in a loop at this line: https://github.com/StochSS/GillesPy2/blob/main/gillespy2/solvers/numpy/ssa_solver.py#L264 However, storing and returning all events that occur is prohibitively expensive and very detrimental to the performance of the solvers.

[Princeofthebow]

Thanks - does the same apply to the C solver?

[briandrawert]

Yes, that same line is here: https://github.com/StochSS/GillesPy2/blob/main/gillespy2/solvers/cpp/c_base/ssa.cpp#L72 It would not be difficult to print out the event times.

[Princeofthebow]

Thank you @briandrawert . If I have time I will try to implement the option. Unfortunately I'm theoretician and tend spend very little time coding but I might do it! Perhaps I could cooperate with @jonstrutz11 ?

[jonstrutz11]

Hi @Princeofthebow, I just ended up using Cayenne which suited my needs.

[dave-doty]

To anyone else finding this thread and looking for a way to simulate discrete chemical reaction networks with all reaction events recorded, we tried Cayenne as suggested by @jonstrutz11. Unfortunately, Cayenne is no longer maintained (and does not work with the latest version of Numpy since it uses the now-missing numpy.int type), but the rebop package does support this: https://github.com/Armavica/rebop/issues/16