jmason42
commented
5 years ago Sure, I have some thoughts.
First, I'll mention that the automated switching between levels of simulation (exact and costly or approximate but cheap) is exactly what tau-leaping does. I personally wouldn't bother trying to unify ODEs and the standard Gillespie SSA because the tau-leaping approach, in the limit that ODEs would be applicable, would still be very fast and give you integral numbers of molecules. If you absolutely need determinism, that's another issue.
Now, as far as ODEs go, translating a stochastic system of reactions into an ODE system is not hard. What's hard is taking the result and mapping it back to integral numbers of molecules. In fact, I bet there are some ugly pathological cases, and even non-pathological cases will still require some careful rounding (and likely even some random sampling). On top of all that, I'm not convinced that these would be any faster to evaluate, given the potential stiffness of things like molecular binding reactions.
I do have some ideas about a quasi-determnistic system that still works on integral numbers of reactions (see #4). In short, I would replace the exponentially distributed events with a degenerate distribution (one in which the event only happens at a prescribed time, effectively the opposite of a "memoryless" distribution). This would be deterministic in all but very specific cases (run conditions in which two exclusive events would occur at the exact same time). It would still need to evaluate one reaction at a time, but usually without any random number generation. However - despite these attractive properties - the physical justification for such a system is weak, particularly in the case of molecules randomly encountering each other.
With all this is mind, you could still consider hybridizing ODEs and exact SSAs. I'd start with the same math as the tau-leaping derivations, but replace the Poisson distributions with degenerate distributions. I'm not exactly sure how far that would take you.
eagmon
It might be useful to switch between stochastic and deterministic simulation based on the time step and molecule counts of each reaction. For example, stochastic complexation is currently finding the same output as its deterministic counterpart. This is because the time step is large, and the counts of most molecules are probably large enough. So at the current state, it would be more efficient to just use a deterministic simulation.
I would think there could be a method for determining whether a reaction should be simulated stochastically or with ODEs, and this can parse out the problem into whatever solver is needed. Ideally, the stochastic method would only be invoked when necessary.
@jmason42, you probably have the best understanding of whether this is feasible, and under what conditions such a switch could be triggered. Do you have thoughts on this?