# Best way to handle (large) chemical reactions
**URL:** https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195
**Category:** New to Julia
**Created:** [October 2, 2017, 12:32pm UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195 "2017-10-02T12:32:36Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [October 2, 2017, 12:32pm UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195/1 "2017-10-02T12:32:36Z")
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](https://github.com/sdwfrost/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
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [October 2, 2017, 2:01pm UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195/2 "2017-10-02T14:01:20Z")
I think you need to say how the chemical reactions are specified.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 2, 2017, 2:49pm UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195/3 "2017-10-02T14:49:09Z")
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](https://github.com/JuliaDiffEq/DiffEqBiological.jl/issues/2)
There’s also BioSimulator.jl to consider:
[https://github.com/alanderos91/BioSimulator.jl](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.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [October 4, 2017, 6:11am UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195/4 "2017-10-04T06:11:31Z")
Thank you for these links, it is very helpful.
---
### Author:  [@erlebach1](https://discourse.julialang.org/u/erlebach1)
#### Post date: [May 26, 2020, 11:47am UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195/5 "2020-05-26T11:47:57Z")
> [@ChrisRackauckas](#):
>
> DiffEqBiological
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.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 26, 2020, 1:01pm UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195/6 "2020-05-26T13:01:05Z")
> [@erlebach1](#):
>
> 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.
We’re getting there: [NetworkSystem · Issue #341 · SciML/ModelingToolkit.jl · GitHub](https://github.com/SciML/ModelingToolkit.jl/issues/341) . I am trying to implement explicit I/O of systems today so that @asinghvi17’s implementation here gives an ODE instead of a DAE.
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [May 26, 2020, 10:52pm UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195/7 "2020-05-26T22:52:30Z")
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…)
```julia
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)
```
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [May 26, 2020, 10:54pm UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195/8 "2020-05-26T22:54:58Z")
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.
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [May 30, 2020, 12:44pm UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195/10 "2020-05-30T12:44:50Z")
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](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 `ReactionSystem`s to ODE/SDE or jump process models:
[https://mtk.sciml.ai/dev/systems/ReactionSystem/](https://mtk.sciml.ai/dev/systems/ReactionSystem/)
---
### Author:  [@erlebach](https://discourse.julialang.org/u/erlebach)
#### Post date: [May 30, 2020, 12:55pm UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195/11 "2020-05-30T12:55:49Z")
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!
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [May 30, 2020, 4:32pm UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195/12 "2020-05-30T16:32:02Z")
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).
---
### Author:  [@erlebach](https://discourse.julialang.org/u/erlebach)
#### Post date: [May 30, 2020, 5:05pm UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195/13 "2020-05-30T17:05:53Z")
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
---
### Author:  [@johnh](https://discourse.julialang.org/u/johnh)
#### Post date: [May 30, 2020, 5:58pm UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195/14 "2020-05-30T17:58:52Z")
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](https://blogs.sciencemag.org/pipeline/archives/2010/02/23/things_i_wont_work_with_dioxygen_difluoride)
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 30, 2020, 6:36pm UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195/15 "2020-05-30T18:36:42Z")
> [@erlebach](#):
>
> 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.
Working on it: [NetworkSystem · Issue #341 · SciML/ModelingToolkit.jl · GitHub](https://github.com/SciML/ModelingToolkit.jl/issues/341)
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [May 30, 2020, 6:46pm UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195/16 "2020-05-30T18:46:49Z")
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.
---
### Author:  [@erlebach](https://discourse.julialang.org/u/erlebach)
#### Post date: [May 30, 2020, 6:55pm UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195/17 "2020-05-30T18:55:36Z")
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.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 30, 2020, 6:58pm UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195/18 "2020-05-30T18:58:32Z")
> [@erlebach](#):
>
> Inference with DiffEqFlux can also be interesting.
Yes, I need to update a few things for that.
---
### Author:  [@ericphanson](https://discourse.julialang.org/u/ericphanson)
#### Post date: [May 30, 2020, 7:07pm UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195/19 "2020-05-30T19:07:27Z")
> [@erlebach](#):
>
> Is it possible to get notices of new messages in this thread even if I don’t participate?
Yes, at the bottom of the thread you can choose “watching” in this dropdown to get notifications:
---
### Author:  [@erlebach](https://discourse.julialang.org/u/erlebach)
#### Post date: [May 30, 2020, 7:16pm UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195/20 "2020-05-30T19:16:36Z")
I meant the thread on GitHub.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 30, 2020, 7:29pm UTC](https://discourse.julialang.org/t/best-way-to-handle-large-chemical-reactions/6195/21 "2020-05-30T19:29:13Z")
Click the subscribe button on the right in the Github issue.