sciml / catalyst.jl Goto Github PK

While direct SBML would be great I suspect it will be a while till someone gets around to that.

As a stop gap we might consider allowing import of simpler formats that can be generated by other tools, such as SBML converters.

In particular, BioNetGen is fairly popular, and has its own language, BNGL, that allows the specification of reaction networks much more compactly than SBML. BioNetGen can also, I believe, import SBML (though I'm not sure what subset of SBML features it currently handles). What's nice about this is they have a greatly simplified "network" file format that is very similar to our DSL, and which they use internally to feed into their solvers. (So I think the BioNetGen workflow is to parse the BNGL/SBML network, expand out the reactions, generate a network file, and feed this file as input into stochastic or ODE solver executables.) The network generation program that generates this reduced file can be run through the commandline. At least an older version of this format is specified in the references to the BioNetGen bible.

So we could consider wrapping the BioNetGen network creation executable with a separate Julia Package (BioNetGen is MIT FWIW), to allow generation of ".net" files, and then add a parser for the ".net" file. The only real challenge I see with the latter is that I think the ".net" files now allow the specification of mathematical expressions for non-standard rate laws. Otherwise the format is quite simple, see the BCR network here. (Note, the generated species names are insanely long as they are coming from a more compact BNGL specification, and encode different complexes / phosphorylation states.)

So I guess the question is

1. Does it make sense to create a package wrapping the BioNetGen network creation executable?
2. Would a ".net" file parser for importing files processed by the wrapper (or standalone BioNetGen) be of interest to anyone else?

Discourse is suddenly overflowing with mention of packages that could be useful for steady-state / bifurcation analysis. We should look into creating interfaces to use them.

PseudoArcLengthContinuation.jl
Bifurcations.jl
HomotopyContinuation.jl

I guess @TorkelE has already started on this to some degree with BRNEquilibrate?

In relation to #1 & my own issues when testing a simple model, I am not sure #7 is entirely correct.

For a true, stochastic, second order reaction with 2 molecules of the same species reacting with each other the reaction rate is
ku[1](u[1]-1) / 2!,
and not
ku[1](u[1]-1).

In the case of 3 molecules the rate would be
ku[1](u[1]-1)*(u[1]-2) / 3!.

This is because the number of ways of pairing 2 molecules of the same species (using binomial coefficients) is
u[1]! / (2! * (u[1]-2)!).
For N molecules this generalizes to
u[1]! / (N! * (u[1]-N)!).
This is a key difference between deterministic and stochastic kinetics.

They apply these corrections in the previously mentioned article, http://stochpy.sourceforge.net/html/userguide_doc.html#stochastic-rate-equations, although without explaining where the constant 0.5 comes from. (Likewise the constant 1/6 in the 3rd order example 3 comes from the factor 1/3!)

There is more information in the original Gillespie papers (e.g. eq 6a in http://www.phys.ufl.edu/~hagen/Readings/gillespie_1977.pdf)

https://github.com/JuliaDiffEq/DiffEqBiological.jl/blob/master/src/reaction_network.jl#L310

That needs to take a save_positions argument. @Gaussia have you been playing with doing this? I know you asked on the Discourse so I assumed you might have this modified already locally?

@isaacsas

We need some way of specifying upregulation and downregulation. We can start by always making it a Hill function of (n,alpha,omega) where alpha A^n / (omega + A^n) multiplies some term. Then we need to allow different forms for specifying when two or more things regulate the same term.

I'm encountering a weird and very specific bug. This occurs when I define a @reaction_network inside a module which itself is precompiled (it took me half a day to pinpoint this one). The definition is ok, but the printout of the type-specification or the SymEngine related fields results in segfaults (and kernel restarts in jupyter).

To reproduce this:

create a file TestBio.jl containing

__precompile__()

module TestBio
using DiffEqBiological
@reaction_network TestReaction begin
    1.0, x --> 0
end
end

In a julia session, enter:

push!(LOAD_PATH, "path/to/your/directory/")
using TestBio

TestBio.TestReaction()

For me, this results in

julia> TestBio.TestReaction()
TestBio.TestReaction(TestBio.#2, Expr[:(+(-x))], SymEngine.Basic[
signal (11): Segmentation fault
while loading no file, in expression starting on line 0
_ZN9SymEngine10StrPrinter5applyERKNS_3RCPIKNS_5BasicEEE at /home/Niklas/programs/anaconda3/envs/conda_jl/lib/libsymengine.so.0.3 (unknown line)
basic_str_julia at /home/Niklas/programs/anaconda3/envs/conda_jl/lib/libsymengine.so.0.3 (unknown line)
toString at /home/Niklas/.julia/v0.6/SymEngine/src/display.jl:4
show at /home/Niklas/.julia/v0.6/SymEngine/src/display.jl:10
print_matrix_repr at ./show.jl:1654
#showarray#263 at ./show.jl:1697
unknown function (ip: 0x7f4339f641fd)
jl_call_fptr_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:339 [inlined]
jl_call_method_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:358 [inlined]
jl_apply_generic at /buildworker/worker/package_linux64/build/src/gf.c:1926
jl_apply at /buildworker/worker/package_linux64/build/src/julia.h:1424 [inlined]
jl_invoke at /buildworker/worker/package_linux64/build/src/gf.c:51
show at ./show.jl:1669
unknown function (ip: 0x7f4339f63c26)
jl_call_fptr_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:339 [inlined]
jl_call_method_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:358 [inlined]
jl_apply_generic at /buildworker/worker/package_linux64/build/src/gf.c:1926
show_default at ./show.jl:140
unknown function (ip: 0x7f4339f56d56)
jl_call_fptr_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:339 [inlined]
jl_call_method_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:358 [inlined]
jl_apply_generic at /buildworker/worker/package_linux64/build/src/gf.c:1926
display at /home/Niklas/.julia/v0.6/OhMyREPL/src/output_prompt_overwrite.jl:8
unknown function (ip: 0x7f4339f56246)
jl_call_fptr_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:339 [inlined]
jl_call_method_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:358 [inlined]
jl_apply_generic at /buildworker/worker/package_linux64/build/src/gf.c:1926
display at ./REPL.jl:125
unknown function (ip: 0x7f4339f55f76)
jl_call_fptr_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:339 [inlined]
jl_call_method_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:358 [inlined]
jl_apply_generic at /buildworker/worker/package_linux64/build/src/gf.c:1926
display at ./multimedia.jl:218
unknown function (ip: 0x7f4339f55e32)
jl_call_fptr_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:339 [inlined]
jl_call_method_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:358 [inlined]
jl_apply_generic at /buildworker/worker/package_linux64/build/src/gf.c:1926
do_call at /buildworker/worker/package_linux64/build/src/interpreter.c:75
eval at /buildworker/worker/package_linux64/build/src/interpreter.c:242
eval_body at /buildworker/worker/package_linux64/build/src/interpreter.c:539
jl_toplevel_eval_body at /buildworker/worker/package_linux64/build/src/interpreter.c:511
jl_toplevel_eval_flex at /buildworker/worker/package_linux64/build/src/toplevel.c:571
jl_toplevel_eval_in at /buildworker/worker/package_linux64/build/src/builtins.c:496
eval at ./boot.jl:235
unknown function (ip: 0x7f435674ed2f)
jl_call_fptr_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:339 [inlined]
jl_call_method_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:358 [inlined]
jl_apply_generic at /buildworker/worker/package_linux64/build/src/gf.c:1926
print_response at ./REPL.jl:144
unknown function (ip: 0x7f4339f55608)
jl_call_fptr_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:339 [inlined]
jl_call_method_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:358 [inlined]
jl_apply_generic at /buildworker/worker/package_linux64/build/src/gf.c:1926
print_response at ./REPL.jl:129
unknown function (ip: 0x7f4339f55008)
jl_call_fptr_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:339 [inlined]
jl_call_method_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:358 [inlined]
jl_apply_generic at /buildworker/worker/package_linux64/build/src/gf.c:1926
do_respond at ./REPL.jl:646
unknown function (ip: 0x7f4339f53771)
jl_call_fptr_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:339 [inlined]
jl_call_method_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:358 [inlined]
jl_apply_generic at /buildworker/worker/package_linux64/build/src/gf.c:1926
do_call at /buildworker/worker/package_linux64/build/src/interpreter.c:75
eval at /buildworker/worker/package_linux64/build/src/interpreter.c:242
eval_body at /buildworker/worker/package_linux64/build/src/interpreter.c:539
jl_toplevel_eval_body at /buildworker/worker/package_linux64/build/src/interpreter.c:511
jl_toplevel_eval_flex at /buildworker/worker/package_linux64/build/src/toplevel.c:571
jl_toplevel_eval_in at /buildworker/worker/package_linux64/build/src/builtins.c:496
eval at ./boot.jl:235
unknown function (ip: 0x7f435674ed2f)
jl_call_fptr_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:339 [inlined]
jl_call_method_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:358 [inlined]
jl_apply_generic at /buildworker/worker/package_linux64/build/src/gf.c:1926
run_interface at ./LineEdit.jl:1583
unknown function (ip: 0x7f43567c67bf)
jl_call_fptr_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:339 [inlined]
jl_call_method_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:358 [inlined]
jl_apply_generic at /buildworker/worker/package_linux64/build/src/gf.c:1926
run_frontend at ./REPL.jl:945
run_repl at ./REPL.jl:180
unknown function (ip: 0x7f4339f17d02)
jl_call_fptr_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:339 [inlined]
jl_call_method_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:358 [inlined]
jl_apply_generic at /buildworker/worker/package_linux64/build/src/gf.c:1926
_start at ./client.jl:413
unknown function (ip: 0x7f43567a0d28)
jl_call_fptr_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:339 [inlined]
jl_call_method_internal at /buildworker/worker/package_linux64/build/src/julia_internal.h:358 [inlined]
jl_apply_generic at /buildworker/worker/package_linux64/build/src/gf.c:1926
unknown function (ip: 0x401bc8)
unknown function (ip: 0x401612)
__libc_start_main at /build/eglibc-SvCtMH/eglibc-2.19/csu/libc-start.c:287
unknown function (ip: 0x4016bc)
Allocations: 7515063 (Pool: 7513257; Big: 1806); GC: 15
[1]    11906 segmentation fault (core dumped)  julia

Note that it is the printout of the f_symfuncs field that causes the error. If you suppress the printout using a semi-colon it all works.

Upon further testing I found that while @ode_def of ParameterizedFunctions.jl does not immediately cause the same error upon instantiation it does error when you print out the symfuncs field.

I have tried this on both Ubuntu and OSX. The problem goes away if you stop precompiling your module.

Does anyone know what witchcraft this is?

Right now a reaction network can be declared as:

@reaction_network(name, ex::Expr, p...)
@reaction_network(name, scale_noise, ex::Expr, p...)
@macro reaction_network(ex::Expr)

with the last one being a special case, returning an error message that the DSL have been updated.

Two things I have been thinking about:

• In some situations a specific type is not needed
• In #19 I suggested updating to affine noise, like (noise_mult, noise_add)
  Now the structure of all the possible inputs: The custom type, the noise, the reaction network, and the parameters are all different. It feels unnecessary to require to input them al.

While right now the only two correct ways to create a reaction network is:

rn = @reaction_network rnType begin
    1.0, X --> Y
end
rn = @reaction_network rnType η begin
    1.0, X --> Y
end

with

rn = @reaction_network begin
    1.0, X --> Y
end

yielding a deprecating message. The parameters are optional. I think it would be an improvement to make all these kinds of input possible:

rn = @reaction_network begin
    1.0, X --> Y
end
rn = @reaction_network rnType begin
    1.0, X --> Y
end
rn = @reaction_network rnType (noise_mult,1.0) begin
    1.0, X --> Y
end
rn = @reaction_network (noise_mult,1.0) rnType begin
    1.0, X --> Y
end
rn = @reaction_network (noise_mult,noise_add) begin
    1.0, X --> Y
end

(the noise may be declared as a parameter or a constant value)

Again with the possibility of adding parameters at the end. It would be easy to make the program check what input is given and put it in the right box. The biggest issue is if we want to enable the declaration of networks with neither noise or type, in which case we would have to remove the deprecation warning.

89d4c82

It was having issues on v0.7 so I disabled it. We can revisit this after the v0.7 upgrade. @Gaussia

Now that @TorkelE has been kind enough to put in Jacobian support I was thinking a bit about how we might want to handle f and the Jacobian for larger systems (where we don't want to make a dense matrix).

It seems like one approach would be the following: let S denote the net stoichiometry matrix and R(u) the vector of rate laws. Then

du = S * R(u)

so we could simply evaluate R(u) as a vector, and then apply the stoichiometry for each reaction to fill in du. This could be done with a sparse matrix to store S, or some other representation. To make this more automatic I was thinking we could just generate a mutable struct that stores S and the cache vector for R(u), and supplies a functor for evaluating the ODE rhs.

Similarly, for the Jacobian we could store a sparse D[R(u)] and multiply it against S to get D[u].

Does this sound reasonable, or are there other approaches that might be more performant for large systems? The approach above does avoid calculating a rate expression more than once per reaction.

Edit: Or perhaps we should just provide a function that gives the action of the Jacobian on a vector and restrict the underlying linear solvers that can be used.

I think it would be interesting to be able to define chemical reaction networks in a more biological fashion. For example:

rs = @reaction_network begin
  A + B --> C
  2C --> A
end

and have that build the reaction set. Combined with #1, we can make this language expressive enough to include other types of terms like Hill functions, and then really extend the ease of use and have it build ODEs, SDEs, and discrete stochastic equations automatically. This also wouldn't be too difficult, it would just require simple parsing to the enventual API like what ParameterizedFunctions.jl does (maybe auto-calculate things like Jacobians etc. for the ODEs/SDEs as well).

I am stock with an unreasonable solution that I get from Julia ode solver. You might be the right person to ask this. Would you please take a look at this problem? I am confused in working with Julia ode solver.

I want to use Julia solver for the chemical reaction below with k=0 and v=1 parameters. which literally that there is only one reaction (v, L1+R --> C0). But when I use the solver, it surprisingly giving an answer with increasing amount of L2.

(k, v), L1 + R <--> C0

k, L2 --> D0
mail.pdf

So there exists methods to efficiently solve certain affine noise SDE problems (http://juliadiffeq.org/2018/04/15/StableSDE.html). This is cool (and useful) and should be included in this. Admittedly this depends a little bit of the state of those solver and the Modelling Toolkit (which seems to have started getting some attention again?).

If integrating this into the modelling toolkit, and getting affine noise problems "for free" seems far off we could enable the @reaction_network macro to generate the corresponding equations. We could have a flag to the construction like:

SDEProblem(model,u0,tspan,p,affine=true)

to designate that we want to create this type of problem instead. The multiplicative and additive noise could either be taken as the two last parameters in the p vector, or given as a tupple:

SDEProblem(model,u0,tspan,p,affine=true,(k,m))

Just closed the dormant PR we have had on this one for a while (#38). Don't fully remember why it was skipped, but there were some issue with it. Opened an issue on the topic until we revisit it.

The tag name "3.5.2" is not of the appropriate SemVer form (vX.Y.Z).
cc: @isaacsas

I was trying out the @reaction_network macro from the manual, and while the first few examples work, the ones with parameters do not. The errors I get are as follows:

julia> Pkg.status("DifferentialEquations")
 - DifferentialEquations         4.1.0

julia> rn = @reaction_network rType begin
           (kB, kD), X + Y ↔ XY
       end kB, kD
ERROR: MethodError: DiffEqBiological.@reaction_network(::Symbol, ::Expr, ::Expr) is ambiguous. Candidates:
  @reaction_network(name::ANY, scale_noise::ANY, ex::Expr, p...) in DiffEqBiological at /Users/me/.julia/v0.6/DiffEqBiological/src/reaction_network.jl:58
  @reaction_network(name::ANY, ex::Expr, p...) in DiffEqBiological at /Users/me/.julia/v0.6/DiffEqBiological/src/reaction_network.jl:53
Possible fix, define
  @reaction_network(::ANY, ::Expr, ::Expr, ::Vararg{Any,N} where N)

julia> rs = @reaction_network begin
         c1, X --> 2X
         c2, X --> 0
         c3, 0 --> X
       end c1 c2 c3
ERROR: MethodError: no method matching @reaction_network(::Expr, ::Symbol, ::Symbol, ::Symbol)
Closest candidates are:
  @reaction_network(::ANY, ::Expr, ::Any...) at /Users/me/.julia/v0.6/DiffEqBiological/src/reaction_network.jl:53
  @reaction_network(::ANY, ::ANY, ::Expr, ::Any...) at /Users/me/.julia/v0.6/DiffEqBiological/src/reaction_network.jl:58
  @reaction_network(::Expr) at /Users/me/.julia/v0.6/DiffEqBiological/src/reaction_network.jl:64

Many SSAs (Gibson-Bruck, partial propensity methods, RSSA...) require more reaction network info than is easily available from the type generated by the @reaction_network macro. It would be great if the generated type could also include fields that provide a mapping from reaction id to the coefficients of reactant species, and a mapping from the reaction id to the coefficients of product species. Such info is needed to determine which reactions need their propensities updated when another reaction occurs or a species changes its population size.

For a bit more detail, many of these other SSAs require knowing things like

Reactants( rx id ) = set of reactants for a given reaction
Products( rx id) = set of products for a given reaction
DependentRxs( rx id) = set of reactions with jump rates depending on species populations that "rx id" modifies.

All these types of sets can be built from the two mappings I mentioned.

Sam

I would often find it useful to visualise a network with an interaction graph. It should be possible to provide a field with sufficient information for such a graph to be automatically generated.

We could either generate the graph according to some pre-existing graph language (LightGraphs.jl?), or just output some fairly raw data and let that processing be done elsewhere (like the expressions and Latexify.jl)

This is a fairly rough idea, and I have not done much research on it, but it would be cool.

If you have a simple reaction System X <--> Y then it requires you to write
k1, X --> Y
k2, Y --> X
if there were some way to rather write
(k1,k2), X<-->Y
that would save a lot of space (especially for models were most reactions goes in both directions). In this case it seems like you cannot use the arrow "<-->" in expressions. However the possibility of reducing these reactions to a single line would be very useful.

This code:

using DifferentialEquations
rn = @reaction_network begin
    p1*p2, X --> Y
end p1 p2
p = [1.,2.]
dProb = DiscreteProblem([10,10],(0,10.),p)
jProb = JumpProblem(dProb,Direct(),rn)

generates the error

UndefVarError: p1 not defined
in top-level scope at base/none
in JumpProblem at N4fWp/src/problem.jl:7
in #JumpProblem#43 at N4fWp/src/problem.jl:18
in network_to_jumpset at N4fWp/src/massaction_jump_utils.jl:84
in make_majump at N4fWp/src/massaction_jump_utils.jl:65
in eval at N4fWp/src/DiffEqBiological.jl:3 
in eval at base/boot.jl:319 
in top-level scope at base/none

It seems to have to do with multiplying p1 and p2 in the reaction rate and creating the jump problem.
e.g. these pieces of code works fine:

using DifferentialEquations
rn = @reaction_network begin
    p1*p2, X --> Y
end p1 p2
p = [1.,2.]
ODEProblem(rn,[10.,10.],(0.,10.),p)

and

rn = @reaction_network begin
    p1, X --> Y
end p1 p2
p = [1.,2.]
dProb = DiscreteProblem([10,10],(0,10.),p)
jProb = JumpProblem(dProb,Direct(),rn)

Should the rates produced by build_jumps_from_reaction be modified for higher order reactions? For example, see the table on page 5 of this paper.

Ref: https://discourse.julialang.org/t/package-compatibility-caps/15301

I am unsure whenever there is an easy solution to this (but it should be biologically justified).

For a reaction with effect on a specific rate r and stoichiometric change s the deterministic contribution is
sr
with a stochastic equation this becomes
sr + ssqrt(r)N(0,1)
however, it should be possible under some circumstances for
sr < |ssqrt(r)*N(0,1)|
which would lead to negative production (or if s is negative, positive degradation).

For larger values (when SDEs should be run) I do not think this should be an issue, but I'd still be interested in thinking about if there is some way to prevent this (without major slowdown, which might be required and would make it not worth it).

For larger networks it might be good to provide functionality to cache dependency graphs, species to jump mappings, etc. These are currently re-generated each time a jump problem is called.

It would be nice to automatically be able to get out the moment equations according to some moment-closure scheme.

https://ifisc.uib-csic.es/users/raul/CURSOS/SP/Introduction_to_master_equations.pdf
http://web.mit.edu/biophysics/sbio/PDFs/stochastics.pdf
https://www-liphy.ujf-grenoble.fr/pagesperso/bahram/Master_Eq/master.pdf
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3176737/pdf/nihms252393.pdf
https://www.omicsonline.org/open-access/chemical-master-equation-empirical-moment-closure-2329-6577-1000155.pdf
http://www.pnas.org/content/110/35/14261
https://www.sciencedirect.com/science/article/pii/S0096300306000117
https://arxiv.org/pdf/1608.06582.pdf

https://arxiv.org/pdf/1509.00612.pdf
https://arxiv.org/pdf/1608.06582.pdf

Might just be making some simple mistake with things I do not understand, but I have been trying to save simulation data from reaction networks using JLD2, and encounter problems:

using DifferentialEquations, JLD2
rn = @reaction_network begin
    (1,10), X ↔ ∅
end
sol = solve(ODEProblem(rn,[20.],(0.,10.)),Rosenbrock23())
JLD2.@save "out.jld2" sol

gives a big error

cannot write a pointer to JLD file

Stacktrace:
 [1] h5fieldtype(::JLD2.JLDFile{JLD2.MmapIO}, ::Type{Ptr{Nothing}}, ::Type, ::Type{Val{true}}) at /home/user.name/.julia/packages/JLD2/KjBIK/src/data.jl:994
 [2] commit_compound(::JLD2.JLDFile{JLD2.MmapIO}, ::Array{Symbol,1}, ::DataType, ::Type) at /home/user.name/.julia/packages/JLD2/KjBIK/src/data.jl:185
 [3] h5type(::JLD2.JLDFile{JLD2.MmapIO}, ::Type, ::SymEngine.Basic) at /home/user.name/.julia/packages/JLD2/KjBIK/src/data.jl:163
 [4] h5type at /home/user.name/.julia/packages/JLD2/KjBIK/src/data.jl:168 [inlined]
 [5] write_dataset at /home/user.name/.julia/packages/JLD2/KjBIK/src/datasets.jl:521 [inlined]
 [6] write_ref_mutable at /home/user.name/.julia/packages/JLD2/KjBIK/src/datasets.jl:526 [inlined]
 [7] write_ref at /home/user.name/.julia/packages/JLD2/KjBIK/src/datasets.jl:534 [inlined]
 [8] h5convert! at /home/user.name/.julia/packages/JLD2/KjBIK/src/data.jl:658 [inlined]
 [9] write_data(::JLD2.MmapIO, ::JLD2.JLDFile{JLD2.MmapIO}, ::Array{SymEngine.Basic,2}, ::Type{JLD2.RelOffset}, ::JLD2.HasReferences, ::JLD2.JLDWriteSession{Dict{UInt64,JLD2.RelOffset}}) at /home/user.name/.julia/packages/JLD2/KjBIK/src/dataio.jl:178

...


 [32] (::getfield(Main, Symbol("##36#37")))(::JLD2.JLDFile{JLD2.MmapIO}) at /home/user.name/.julia/packages/JLD2/KjBIK/src/loadsave.jl:50
 [33] #jldopen#31(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::getfield(Main, Symbol("##36#37")), ::String, ::Vararg{String,N} where N) at /home/user.name/.julia/packages/JLD2/KjBIK/src/loadsave.jl:4
 [34] jldopen(::Function, ::String, ::String) at /home/user.name/.julia/packages/JLD2/KjBIK/src/loadsave.jl:2
 [35] top-level scope at /home/user.name/.julia/packages/JLD2/KjBIK/src/loadsave.jl:48
 [36] top-level scope at In[1]:6

To be honest I have little clue what is going on here and maybe it is an issue for the people over at JLD2, I will have a check with them as well.

Right now we create a set of ConstantRateJumps, possibly including a few VariableRateJumps.

However recently RegularRateJumps have been added. We could add this to the reaction_network DSL as well. Looking at the syntax the difference is that for a RegularRateJumps you put everything into only one function each for rate and affect, it should be no problem doing the same in reaction_network.

We would have to decide which type of JumpProblem should be generated by default though.

In an attempt to handle all input parameters the anonymous function expressions

(du,u,p,t) -> ...

were replaced with

(internal_var___du,internal_var___u,internal_var___p,internal_var___t) -> ...

As a result of this the function do not understand that you mean the time variable if you input t in the reaction rate (as in something time dependent, e.g.

rn = @reaction_network rType begin
    step(t,10,10), X --> Y
end

You can still do this, but then you have to actually write

rn = @reaction_network rType begin
    step(internal_var___t,10,10), X --> Y
end

The simplest solution is to simply make t a forbidden variable again. Alternatively one can check that t have not been gives as a reactant/parameter. If that is the case all instances of t will be replaced by internal_var___t.

The second alternative will not be that complicated. Unless we think it might trick people to think they are using the time t while they are not.

The tag name "v.3.4.0" is not of the appropriate SemVer form (vX.Y.Z).
cc: @ChrisRackauckas

I thought it would be an interesting to start building in parallel a new DSL for biological reactions on top of SciCompDSL.jl. I have some ideas about it an a draft of the data structures/operators that could be used, but I would like to throw it around first.

First, the system I have in mind so far would have the following features:

• Each reaction must have a rate equation which is just an arithmetic expression like the ones supported in SciCompDSL. This rate equation must explicitly use the reactants in the reaction. This may be cumbersome for simple mass action models but it allows to be more general and even go beyond mass action (e.g. Goldman–Hodgkin–Katz flux equation for ion transport)
• There is no need for directionality in the arrows to specify the reaction as the rate equation already includes the logic. So one arrow suffices for all types of reactions.
• Reactants can be defined to be inside compartments as well as reactions. This allows to simulate transport across compartments.
• Reactions may be mixed with explicit differential equations, intermediate variables, etc (i.e. what is already available at SciCompDSL). The ReactionDSL will simple lower the reactions into intermediate expressions and differential equations and leave the rest untouched.
• A reactant may also appear in an explicit differential equations, intermediate variable, etc and the lowering of reactions will just be merged with what is already there.

Some types are needed:

• Compartment: This can be lower to a Parameter or a DependentVariable based on whether it is constant or changes during simulation (e.g. due to cell growth). A compartment does not have to be a volume, it could be a surface too. The role of compartment is to perform transformations when a reactant and reaction are not defined on the same compartment (see above).
• Species: This means "chemical species", not the biological concept. This is lowered to a DependentVariable. Its differential equation is generated by adding up all the rate equations of the reactions it is involved, modified by stoichiometry and compartment (when different).
• Reactant: The basic component of a reaction. This is basically a Species with an stochiometry. In most cases the stoichiometry would be a constant integer, but non-integer stoichiometries exist (e.g. check ATP synthase) and it may also be an idea to make it a parameter for people designing biological systems that don't exist yet (e.g. synthetic biology).
• Reaction is just a list of reactants (organized as substrates and products) and a compartment.

Operators can be overloaded to generate the right types. For example:

(2p + 2q → 2r => C) == k*(2p*2q - 2r/Keq)

The → can be used to generate the Reaction by analyzing the AST generated on each side (i.e. Species would behave like Variables in SciCompDSL).
=> may be used to indicate the compartment of the reaction (by default the Reaction would get the compartment of the first Reactant so this is not always needed). I would really like to use in instead but then we would need to use parenthesis due to operator precedence.
The == would generate an Operation as in SciCompDSL, with a Reaction as first argument. Because of operator precedence, the parenthesis are required. A function would then process all the Operations that have reactions and lower them. The rest of Operation in the model will be ignored. In the end, the lowered and ignored Operations are passed to SciCompDSL to further processing into an ODE, SDE, etc...

I would love if ReactionDSL would eventually have all the features that DiffEqBiological has already and to make the transition as painless as possible, hence my posting of this message here. In particular, I am not sure whether the "single arrow approach + explicit use of reactants in rate equation" will be seen as step backwards for most people (though, again, not all biological models follow mass action). Also, I have only built ODEs in the past, so I am not sure how this would address the other types of models.

@korsbo @Gaussia

It would be nice to implement various algorithms for qualitative analysis of chemical reaction networks. So, this issue is for brainstorming:

1. What kind of qualitative analyses would you like to see? Preferably with a link to papers, if anyone has published it before.
2. What kind of qualitative analyses could you contribute?
3. What kind of support from the existing framework do we need?

I would plan to implement BF16 (other potential references: OM16, OTM17); in short, this takes a network with qualitative kinetics (we only know stoichiometry, and know on which concentrations each reaction rate depends), and outputs a graph that describes which steady state concentration and steady-state reaction rates depend on which parameters of the system, and also generates some bounds on the kind of possible bifurcations from the steady state. This graph is pretty structured: It has a transitivity property that allow nice visualizations.

I would further plan to implement some variant of "Feinberg Shinar concordance" analysis. I recommend using this as a google term; I know no readable paper on this topic that I can recommend (instead of subjecting any of you to the pain of trying to really understand more than the abstract of this, I'd rather recommend waiting until I post some code).

What kind of things do I need from the network:

I would need to extract, for each reaction, stoichiometry and dependencies (on which concentrations does the reaction rate depend).

Concordance analysis would also require monotonicity information (each reaction rate can depend each concentrations either positively, negatively, with unspecified sign or not at all). Unless we do some extensive interval math, I'd guess that this information would need to be user supplied.

For convenience, one would need a nice way to "clamp" concentrations (remove them from the state-space by considering them constant).

If we want to do some more exotic Feinberg-style analyses which only apply to mass-action kinetics, we would also need to know which reactions are supposed to have strict mass action kinetics. This would need to be a user-supplied parameter: Most of the time, you consider mass-action as an approximation of some unfathomably complicated kinetics; sometimes (rarely) you want to consider it physical/mathematical truth.

If we want the comparison between different variants of networks from section 5 of https://arxiv.org/abs/1606.00279, then we would need some some support from the DSL. Not sure how this should look like.

Afterwards, we would like to visualize the output of the analysis, and store it in some usable data structure. We would further like to use it for the steady state solve, and if parameter-tracking / bifurcation analysis ever becomes a thing in diffeq, then we would like to use it there as well.

For the steady state solver, the result looks like the following: We factorize our big, coupled nonlinear system into a block tridiagonal structure; that is, we have a list of lower-dimensional systems f_i(x_1,..., x_k)==0 (each f_i and x_j is multi-dimensional), and a transitive relation isless, such that f_i(x)=f_i({x_k: isless(k,i)}). Obviously this severely constrains the kind of bifurcations or multiple steady states we can have.

I'm trying to translate a large network (~25,000 reactions) from a simple text file into DiffEqBiological. Right now this seems to be too much to handle. For reference, parsing the network, which is in a simple file with lines of the form:

A + B -> C, 1.0

is super quick. I then translate each line to the appropriate form for use in the @reaction_network macro, and construct a full instance of the @reaction_network macro as a string. Then I call eval( Meta.parse(rnstr) ) on this string. The latter seems to take a very large amount of time, and use up a lot of memory. (Up to 17 minutes and 11GB at the moment.)

Does anyone have a suggestion for a better way to do this?

From the documentation examples I tried:

using DiffEqBiological
rn = @reaction_network rType begin
     (kB, kD), X + Y ↔ XY
end kB, kD

and get

ERROR: MethodError: DiffEqBiological.@reaction_network(::Symbol, ::Expr, ::Expr) is ambiguous. Candidates:
  @reaction_network(name::ANY, scale_noise::ANY, ex::Expr, p...) in DiffEqBiological at /Users/isaacsas/.julia/v0.6/DiffEqBiological/src/reaction_network.jl:58
  @reaction_network(name::ANY, ex::Expr, p...) in DiffEqBiological at /Users/isaacsas/.julia/v0.6/DiffEqBiological/src/reaction_network.jl:53
Possible fix, define
  @reaction_network(::ANY, ::Expr, ::Expr, ::Vararg{Any,N} where N)

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.

With ParameterizedFunctions.jl, the variable names are fed to the legends generated by Plots.jl, it would be nice if DiffEqBiological could do the same.

@isaacsas had an interesting point that since concentrations and amounts should always be positive, why not add PositiveDomain on them by default? @TorkelE you had some example problems that were hitting 0: what happens if we try that on your problems?

#87 (comment)

After I updated my packages, the jump equations from @reaction_networks have acquired a bug. I don't know if it is from here or from DiffEqJump.

The issue occurs when I use a function for the reaction rates

using DifferentialEquations
using Plots


model = @reaction_network test1 begin
    (p_x, d_x), 0 ↔ x
    (hill(x, p_y, 1, 1), d_y), 0 ↔ y
    end p_x d_x p_y d_y

p = [100., 1.0, 50., 1.0]
u0 = [0, 0]
tspan = (0., 10)
prob = ODEProblem(model, float.(u0), tspan, p)
odesol = solve(prob)

discrete_prob = DiscreteProblem(u0,tspan,p)
jump_prob = JumpProblem(discrete_prob,Direct(),model)
sol = solve(jump_prob,FunctionMap())


plot(odesol)
plot!(sol)

Here, the y variable remains at 0 with the Jump equations (which it should not do and which it did not use to do).

I could try to pinpoint where the issue was introduced, but I figured that @isaacsas or @ChrisRackauckas may already know what has gone wrong.

I noticed DiffEqBiological.jl supports ODEProblem SDEProblem, and JumpProblem. Are there plans to support SteadyStateProblem as well? Often in biological models built using ODEs, people are interested in solving for the steady-state fluxes, concentrations, etc. after applying some perturbation (e.g. an enzyme overexpression).

I was reading the documentations today after the latest update and noticed that g could be a sparse matrix like:

A = zeros(2,4)
A[1,1] = 1
A[1,4] = 1
A[2,4] = 1
sparse(A)

# Make `g` write the sparse matrix values
function g(du,u,p,t)
  du[1,1] = 0.3u[1]
  du[1,4] = 0.12u[2]
  du[2,4] = 1.8u[2]
end

# Make `g` use the sparse matrix
prob = SDEProblem(f,g,ones(2),(0.0,1.0),noise_rate_prototype=A)

Right now we make a g including all combinations du[i,j] = , which might be quite a few. It should be very possible to instead make this g a sparse matrix with much less entries (e.g. if we have n components being degraded that will add n squared entries to g, only of which n is non zero).

Could this be an advantage? Would there be a non trivial performance gain by using such a matrix. If that is not the case we probably shouldn't bother.

Wouldn't it be nice if we could provide fields with symbolically calculated jacobians, etc., like in ParameretizedFunctions?

Me and @Gaussia are already working on implementing this, but we are unsure of how to best realise the idea. It should be pretty straightforward to actually get the things calculated (with heavy inspiration from ParameterizedFunctions.jl) but I'm wondering how to best present the available expressions to the problem/solver.

The Reaction type that we are working on has a field for a function with the ODE. Let's say that we have an instance of the Reaction type called r. The ODE will be available as r.f. Would it be best to place the symbolically calculated stuff in separate fields of r, or should we just make r.f a subtype of AbstractParameterizedFunctions and place all the stuff in there (e.g. r.f.symjac).

And if not, what do we need to overload in order for the DEProblem/solver to know that the jacobian, etc., exists?

@isaacsas and @ChrisRackauckas Moving some stuff over to a separate conversation.

I plan to update the ReactionStruct to something like this:

struct ReactionStruct
    substrates::Vector{ReactantStruct}
    products::Vector{ReactantStruct}
    dependants::Vector{Symbol}
    rate_org                 #(an expression, but might become just a symbol or number)
    rate_DE                 #(an expression, but might become just a symbol or number)
    rate_SSA               #(an expression, but might become just a symbol or number)
    is_pure_mass_action::Boolean

I want to have as few fields as possible while including everything that is needed.

• rate_org is the rate as inputed in the DSL. This could e.g. be 3.5. The rate needs to be saved in its original form since MassActionJumps needs it as input. T
• rate_DE the rate as used in ODEs and SDEs, modified after mass action (unless ⟾ is used...). Something like 3.5*X^2/2!.
• rate_SSA yet a third form is needed for rates for normal jumps (like RegularJumps and ConstantRateJumps). This is different from the rate_DE since it might modify copy number like 3.5*X*(X-1)/2!.
• dependants is a list of dependants, as required to make dependency graphs.
• is_pure_mass_action is true if a MassActionJump is used.

Unless we need something else, or if there is a way to make some of this redundant this is probably what we will go with. A primary purpose of this is so that @isaacsas knows what he will have to work with for making dependency graphs (you start on that end).

Coding the repressilator:

repressilator = @reaction_network begin
    hillr(P₃,α,K,n), ∅ --> m₁
    (δ,γ), m₁ ↔ ∅
    hillr(P₁,α,K,n), ∅ --> m₂
    (δ,γ), m₂ ↔ ∅
    hillr(P₂,α,K,n), ∅ --> m₃
    (δ,γ), m₃ ↔ ∅
    β, m₁ --> m₁ + P₁
    μ, P₁ --> ∅
    β, m₂ --> m₂ + P₂
    μ, P₂ --> ∅
    β, m₃ --> m₃ + P₃
    μ, P₃ --> ∅
end α K n δ γ β μ

gives:

julia> repressilator.f_func
6-element Array{Expr,1}:
 :((α * K ^ n) / (K ^ n + P₃ ^ n) + -1 * δ * m₁ + γ + 0)
 :((α * K ^ n) / (K ^ n + P₁ ^ n) + -1 * δ * m₂ + γ + 0)
 :((α * K ^ n) / (K ^ n + P₂ ^ n) + -1 * δ * m₃ + γ + 0)
 :(β * m₁ + -1 * μ * P₁)                                
 :(β * m₂ + -1 * μ * P₂)                                
 :(β * m₃ + -1 * μ * P₃)     

While I assume the extra +0 gets optimized away when the actual ODE function is compiled, it would be nice to get it out of the expression since it appears in other contexts (such as displaying the ODEs with Latexify).

This is for the new DiffEq parameter update. I am not sure how to make those tie into the current setup's rate parameters.

I think almost all the remaining Anys should be replaceable by Union{Expr,Symbol,Number}.

DiffEqBiological.jl makes it easy to go from a reaction network (defined using the DSL) to an ODEProblem.

However, large biological models (100s of reactions +) are often encoded in a more convenient format such as SBML. That said, it would be great if we could input a stoichiometric matrix directly instead of listing the reactions in the DSL, as an alternative input. We would also need to input a vector of rates as well.

This would allow those with large biological models to go from a stoichiometric matrix (which can be generated using COBRA tools) directly to an ODEProblem.

I noticed that currently, the reaction networks DSL only takes stoichiometry into account. So it treats A -> 0 and A + A -> A as the same system.

However, because of the law of mass action, the former is just exponential decay (a decay rate proportional to the concentration) with a solution of the form [A_0] e^-kt , while the second is a rational decay (the decay rate is proportional to the concentration squared) with a solution of the form 1 / (k t + [A_0]^-1 ) where k is the reaction rate and A_0 the initial concentration.

The two behave rather differently even from a qualitative point of view. In a system of equations, A would stick around and influence the system for much longer in the second case. If we add a second row A + B -> A at a much smaller rate than the decay rate of A, then in the first case the concentration of B will be negligibly affected, while in the second case it will be significantly affected as the integral of A over time diverges.

Not sure if these are something to be concerned about, but thought I'd report them in any case.

sir_model = @reaction_network SIR begin
           c1, s + i --> 2i
           c2, i --> r
       end c1 c2
p = (0.1/1000,0.01)
u0 = [999.,1.,0.]
tspan = (0.,250.)

julia> @code_warntype ODEProblem(sir_model,u0,tspan,p)
Body::Any
34 1 ─ %1 = invoke DiffEqBase.:(#ODEProblem#302)($(QuoteNode(Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}()))::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, _1::Type, _2::Function, _3::Array{Float64,1}, _4::Tuple{Float64,Float64}, _5::Tuple{Float64,Float64})::Any
   └──      return %1         

solve for ODEs seems to be fully concretely typed. For SDEProblem both problem construction and solve give Anys:

julia> @code_warntype SDEProblem(sir_model,u0,tspan,p)
Body::Any
2 1 ─ %1 = (getfield)(args, 1)::Tuple{Float64,Float64}                   │      
  │   %2 = (getfield)(args, 2)::Tuple{Float64,Float64}                   │      
  │   %3 = (Base.getfield)(rn, :g)::Function                             │╻      #SDEProblem#44
  │   %4 = (Base.getfield)(rn, :p_matrix)::Array{Float64,2}              ││┃      getproperty
  │   %5 = %new(NamedTuple{(:noise_rate_prototype,),Tuple{Array{Float64,2}}}, %4)::NamedTuple{(:noise_rate_prototype,),Tuple{Array{Float64,2}}}
  │   %6 = DiffEqBiological.SDEProblem::Type{SDEProblem}                 ││     
  │   %7 = invoke getfield(Core, Symbol("#kw#Type"))()(%5::NamedTuple{(:noise_rate_prototype,),Tuple{Array{Float64,2}}}, %6::Type{SDEProblem}, _2::Function, %3::Function, _3::Array{Float64,1}, %1::Tuple{Float64,Float64}, %2::Tuple{Float64,Float64})::Any
  └──      return %7                                                     │      
julia> @code_warntype solve(prob,EM(), dt=.01)
Body::Any
 1 ─ %1 = %new(Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol},NamedTuple{(:dt,),Tuple{Float64}}}, #temp#, (:dt,))::Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol},NamedTuple{(:dt,),Tuple{Float64}}}
 │   %2 = (getfield)(args, 1)::Core.Compiler.Const(EM{true}(), false)                 │  
 │   %3 = invoke DiffEqBase.:(#solve#442)(%1::Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol},NamedTuple{(:dt,),Tuple{Float64}}}, _3::Function, _4::SDEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Float64,Float64},Nothing,SDEFunction{true,SIR,getfield(Main, Symbol("##9#25")),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Array{Symbol,1},Nothing},getfield(Main, Symbol("##9#25")),Nothing,Array{Float64,2}}, %2::EM{true})::Any
 └──      return %3     

for jumps only problem construction gives Any

julia> @code_warntype JumpProblem(dprob,Direct(),sir_model)
Body::Any
7 1 ─ %1 = invoke DiffEqBiological.:(#JumpProblem#45)($(QuoteNode(Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}()))::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, _1::Type, _2::DiscreteProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Float64,Float64},DiscreteFunction{true,getfield(DiffEqBase, Symbol("##270#271")),Nothing},Nothing}, _3::Direct, _4::SIR)::Any
  └──      return %1