Open sophiejwalton opened 4 years ago
This is in general a difficult problem to solve, but it can be solved if it is actually an equilibrium (as opposed to more generically a steady state). In the equilibrium case, I have written a package to do that, which will be released soon. The devs here and @sophiejwalton know where to find me if you want to check it out.
@sophiejwalton : Are you referring to computations of concentrations that are at quasi-steady state (as discussed here - Lecture 2 BE 240 / or if a species is at equilibrium ? Can you give an example.
@justinbois Can you give a link to the package you are referring to? I would love to check it out. I have been writing a package to do QSSA based model reductions in an automated manner - it's available here - in a very early stage.
@justinbois this is a very cool feature we were hoping to include in Bioscrape at some point - what is your package written in? The thought we have is to have "Algebraic Rules" for equilibrium / quasi-steady-state which can be solved. For complex balanced CRNs in the deterministic regime, we know there is a well defined solution. In other cases, there may be multiple solutions which presents a complication. Also we aren't sure how to extend to the stochastic case where we want to sample from the equilibrium distribution somehow...
The package is written in Python, though there is legacy C code for it as well. The Python code is almost entirely JITted using Numba, so the first execution is slow, but every execution thereafter is fast.
It solves the problem of coupled equilibria: Given a set of reversible chemical reactions, their equilibrium constants, and the initial concentrations of all chemical species, compute the equilibrium concentration of all species assuming the law of mass action holds.
The algorithm is provably globally convergent, and it's quite fast.
We plan to write it up as a paper and release it soon.
I look forward to seeing the paper! I'd love to try to convince you to include this as cython code integrated with bioscrape.
Do you have any sense if there is a stochastic analog to this algorithm in other words a way to quickly sample the equilibrium distribution?
EQTK is now released if you want to use it. https://eqtk.github.io/.
Thanks @justinbois - conceptually, this looks really cool! Can you direct me to anything more theoretical about how it works?
The original paper that lays out the principles is here: https://doi.org/10.1137/060651100
We are working on another paper now. I can forward you via PDF further notes.
This is Sophie from BE 240. It would be nice to have a feature to compute exact concentrations of complexes at equilibrium for complex binding functions.
Like what if A+ B => C C + C => D or Please let me know if this makes sense.