tahorst
commented
6 years ago With implementing the reverse reactions with Gillespie, do you think that would increase the number of reactions we would have to model with this change? My understanding is it essentially goes to completion now so there isn't really any reverse but if we do allow for those stochastic reactions to happen in reverse we might have more than 2 molecules that we're solving for in the later time steps. I would guess it wouldn't be too much more so you might be right that there would be minimal iterations of Gillespie to do but still something good to check on. It does sound like a good idea for both speed and modeling accuracy to implement this though.
jmason42
In looking at the ODEs used ultimately for the equilibrium and two component system processes, I noticed that these ODEs are stiff and effectively impossible to solve without a specialized solver. This is because the modeled rates of reaction e.g.
A + B <-> Care predominantly forward; that is, the rate of the reactionA + B -> Cis much higher than the rate of reactionA + B <- C. This creates the numerical problems that make this problem difficult to solve.My proposal is to redesign this problem as a stochastic simulation - either exact (Gillespie) or approximate (tau-leaping). There are a couple reasons for this; for one, the number of molecules and reactions is very small. In the first time-step, for example, the two-component system process only models the reactions between about a thousand molecules. In later time steps, it's only about 2 molecules. That's two reactions to solve for in a given time-step, meaning two iterates of the Gillespie algorithm. Very cheap compared to solving a system of stiff ODEs, and then rounding the result to make the numbers into integers. Additionally, this would allow for rare, stochastic un-binding events. Right now the state of the system is deterministic.
I personally don't have time to dedicate to this, but I would be happy to assist anyone who is interested.