Hi,
I am wondering what is the best way to handle chemical reactions for SSA simulations. For example, given a set of chemical reactions, I want to extract the transition matrix and rate function in order to use Gillespie.jl. One could use the python package of Stochpy or pyxlr8, but I was thinking about a “julian” compact way of doing this.
Thank you for your help / sharing your code,
Best regards
I think you need to say how the chemical reactions are specified.
Well you probably don’t want to use a transition matrix for a large sparse set of reactions. DifferentialEquations does it all using a sparse representation so that should do better on large networks for a small cost on smaller networks.
There is also an undocumented DSL for handling this.
https://github.com/JuliaDiffEq/DiffEqBiological.jl/issues/2
There’s also BioSimulator.jl to consider:
https://github.com/alanderos91/BioSimulator.jl
I do plan on wrapping these all onto the DiffEq API if they keep developing. If there’s something else you’re looking for like some file parsing, feel free to suggest in DiffEqBiological.jl. I am unaware of a common format and if we should support it.
Thank you for these links, it is very helpful.
So far, these equations are all defined at a single point. Is there software that considers network connectivity? For example, a graph (G,E) with equations defined on the edges. Each edge would have its set of reaction equations. The reactants would live at the notes.
Thanks.
We’re getting there: NetworkSystem · Issue #341 · SciML/ModelingToolkit.jl · GitHub . I am trying to implement explicit I/O of systems today so that @asinghvi17’s implementation here gives an ODE instead of a DAE.
Not sure what you are looking for, but here is how you can create a ModelingToolkit ReactionNetwork easily from a LightGraphs DiGraph, convert it to a JumpSystem and then simulate it. (I wrote this pretty quick, so hopefully there are no bugs, but it should give the idea…)
using ModelingToolkit, DiffEqBase, DiffEqJump, LightGraphs
nedges = 7
nverts = 4
dg = DiGraph(nverts, nedges)
@parameters t k[1:nedges]
@variables u[1:nverts](t)
rxs = [Reaction(k[i],[u[src(e)]], [u[dst(e)]]) for (i,e) ∈ enumerate(edges(dg))]
rs = ReactionSystem(rxs, t, u, k)
js = convert(JumpSystem, rs)
# each vertex's value is a random integer in 1...100
u0map = [u[i] => rand(1:100) for i ∈ eachindex(u)]
pmap = [k[i] => rand() for i ∈ eachindex(k)]
tspan = (0.0,1.0)
dprob = DiscreteProblem(js, u0map, tspan, pmap)
jprob = JumpProblem(js, dprob, NRM())
sol = solve(jprob, SSAStepper())
using Plots
plot(sol)
And just a note, but for a system that small Direct() is probably the Gillespie method to use in DiffEqJump. For large, sparse networks there are a bunch of other methods that can perform better there.
I see you deleted your post, but just in case you still have those questions:
NRM is one of the DiffEqJump exact SSAs (i.e. Gillespie method variants). Direct is the classical Direct method of Gillespie. You can see a bit of information about the (exact) methods we’ve implemented at:
https://docs.sciml.ai/dev/types/jump_types/#Constant-Rate-Jump-Aggregators-1
For small systems generally Direct is the fastest. For large systems one of SortingDirect, DirectCR or RSSA is usually fastest, depending on the network topology.
Reaction is defined in ModelingToolkit, where we recently added a ReactionSystem type one can define, and methods to convert ReactionSystems to ODE/SDE or jump process models:
Thank you! I deleted the account once I figured out the answer. Perhaps it was bad form. Your reply provides details I did not know. One lingering question is this. Is there any need for DiffEqBiological.jl? The Reactions routine seems tomdo the job.
Finally, I am implementing an SIR model and have equations for S, I, and R on each edge. I assume that all I have to do is use Reaction three times on each edge?
Thanks!
DiffEqBiological is currently being revamped to serve two purposes. It will provide an easier interface for specifying ModelingToolkit ReactionSystems, avoiding having to define each Reaction directly. It will also provide an API for querying and modifying the generated ReactionSystem (i.e. features like merging two togther and such). If you are comfortable just building the ReactionSystem directly there won’t be any real functionality that is lost from bypassing the reaction_networik macro.
For a graph-based SIR model I’d do exactly what you mention. Just define array variables for the S, I and R species, each index corresponding to the value of a species at a vertex of the graph. Then define three reactions per edge to represent migration of each species between vertices (i.e. S_i \to S_j). You can also define reactions for each species at the same vertex to have the normal infection dynamics conversion between S_i, I_i and R_i locally at vertex i (i.e. S_i+I_i \to 2I_i and such).
Thank you, Sam! I will try this out (SIR model)
It would be nice though if one could define an equation set once and have a method that applies it to all edges. Ideally, apply an arbitrary function with arbitrary arguments to all the edges of a graph, or to an edge set. That would further simplify construction.
Gordon
Best way to handle large chemical reactions? Run away! Run away!
Satan’s kimchi
https://blogs.sciencemag.org/pipeline/archives/2010/02/23/things_i_wont_work_with_dioxygen_difluoride
Yeah, such functionality would be great, and for my own research would be helpful for defining spatial Gillespie models. It is in the pipeline, but nothing concrete has been setup to enable it yet.
Is it possible to get notices of new messages in this thread even if I don’t participate?
Inference with DiffEqFlux can also be interesting.
Ever heard of Multilayer networks? Any support you know of?
I do not need it now.
Yes, I need to update a few things for that.
I meant the thread on GitHub.
Click the subscribe button on the right in the Github issue.
