Build ODEs and SDEs from reactions

[ChrisRackauckas]

There should be overloads for ODEProblem(rs::ReactionSet,u0,tspan) and SDEProblem(rs::ReactionSet,u0,tspan) which allows the user to define ODEs and SDEs directly from a reaction set which builds the mass action model.

[gajomi]

I think this makes a lot of sense and would make quite a few people happy. One key issue would be to make sure that one has all the information needed to deduce the appropriate ODE/SDE. By "information" I refer to not only data int the ReactionSet, but also information about the scale (so that one can convert numbers into densities and visa versa) as well as the approximation scheme itself, which might yield different continuous rate functions for the same reaction model. For example, from a user perspective, one might like to make a basic distinction between a system with effectively point like reacting particles versus a system with crowding (e.g. on the surface of a cell or along a polymer) where one one has "null" or empty particles.

The above might be more of a rant than a useful contribution, but one concrete thign that would also help would be to complie a list of example reaction set/sde/ode triplets, which might help flesh out what extra information may or may not be needed to make this happen. To get the ball rolling I would offer the examples in this review paper: https://arxiv.org/pdf/q-bio/0501022.pdf.

[ChrisRackauckas]

I was just going to have it use the convergence theorems from Gillespie and create the large N ODE and SDE approximations (derived via the Fokker-Plank equation). The next step I was going to take after that was interactions and regulation, allowing Hill functions and specifying how they combine. I am not sure how things like crowding fit into this, though there would be no issues with having a separate DSL if necessary.

[gajomi]

Cool. I think I (sort of) follow.

But my (potentially incorrect) understanding is that there isn't just one way to go from the reaction rates to the ODE/SDE, but rather multiple methodologies that may yield different results (e.g. pair expansions, Van Kampen approximation, etc.). So I guess I am trying to say that at minimum there should be some kind of documentation as to the methodology employed here.

As for my remark about "crowding" (probably not the right word), I think the Hill function feature you mentions interacts directly with this and the aforementioned methodological choices. In terms of the microscopic dynamics, the Hill function from both (a) having a finite number of binding sites available on a surface/polymer/whatever as well as (b) an argument that the binding/unbinding process is fast relative to whatever else is happening (dimerization, degradation, etc.). These are interesting together in the context of your proposal because depending on how one considers the associated stochastic fluctuations one can end up with different noise terms in the Fokker Plank approximation.

Anyway, sorry for rantlike wall of text. I hope it doesn't come off as a shotgun criticism. But I am keenly interested in this sort of thing and may be interested in contributing if I could get a better sense for the scope of the problems you are interested in solving with it.

[ChrisRackauckas]

No worries at all. What I am alluding to is that there is a pretty standard way that phenomenological models of gene regulatory networks are developed with mass action laws which just involves mass action and simple Hill function regulation. I want to get that first before considering more.

[TorkelE]

Hello, Since I am spending most of my time modelling reaction networks I have had my own program for reading reaction networks. It is quite similar to yours but with some additional function. The syntax is something like:

network = @read_network begin
    1.0, X --> Y
end

were network is a structure like:

struct network_model
    f::Function
    g::Function
    jumps::Tuple{ConstantRateJump,Vararg{ConstantRateJump}}
    p_matrix::Array{Float64,2}
end

Next I can use that to solve ODE/SDE or make Guilespie simulations using

prob_det = ODEProblem(network.f,u0_New,tspan)
prob_stoch = SDEProblem(network.f,network.g,u0,tspan,noise_rate_prototype=network.p_matrix)
prob_jump = JumpProblem(prob_disc,Direct(),network.jumps...)

I have developed it for my own use, but if you are interested I can share the code with you, there might be something useful for you. The module itself is just under 200 lines and also contain some additional functionality: Reactions in both directions: (1.0,2.0), X ⟷ Y Constants or functions in the reaction rates: (hill(Y,2,3.0,1.5),kD), X ⟷ Y You can also use an open arrow to make it ignore mass kinetics and simply use whatever reaction rate you have designated 2.0, X ⇒ Y will have the reaction rate 2.0, als opposed to 2.0[X]. I've also added support for all types of arrows I could find.

I am a bit unsure what you are looking for, or if you have plans to expand this part of the package. But if you do I am happy to share what I got, certainly some of it might be useful.

[ChrisRackauckas]

Wow, that sounds perfect. That's what I've been looking to find the time to implement:

Constants or functions in the reaction rates:

That would solve https://github.com/JuliaDiffEq/DiffEqBiological.jl/issues/13

Reactions in both directions:

That would solve https://github.com/JuliaDiffEq/DiffEqBiological.jl/issues/11

functions in the reaction rates:

That would solve https://github.com/JuliaDiffEq/DiffEqBiological.jl/issues/8

I would greatly appreciate a PR! This sounds like the solution to our problems. Only thing I would change is:

prob_stoch = SDEProblem(network.f,network.g,u0,tspan,noise_rate_prototype=network.p_matrix)

that should just get a dispatch

prob_stoch = SDEProblem(network,u0,tspan)

Please PR 👍

[TorkelE]

Will do, cheers :) Just make some adaptations to fit the code in your framework. Also my friend, who unlike me have some experience in Julia, is giving my some tips about Julia standard. Will PR soon.

[ChrisRackauckas]

No worries about making it perfect. I can help massage the code as well during review.

[ChrisRackauckas]

Fixed by https://github.com/JuliaDiffEq/DiffEqBiological.jl/pull/15