SciML / ReactionNetworkImporters.jl

Julia Catalyst.jl importers for various reaction network file formats like BioNetGen and stoichiometry matrices
https://docs.sciml.ai/ReactionNetworkImporters/stable/
MIT License
26 stars 9 forks source link

support BioNetGen function blocks and more inline function types #4

Open isaacsas opened 5 years ago

isaacsas commented 5 years ago

Currently function blocks are not supported. This would be great to add but note, this requires allowing observables within the functions.

pascalschulthess commented 2 years ago

Yes. functions and groups would be great!

isaacsas commented 2 years ago

Do you have an example that doesn’t currently work?

pascalschulthess commented 2 years ago

I'm currently looking into this: https://github.com/labsyspharm/marm1-supplement/blob/master/resources/MARM1.bngl which may be way to complex but it does have both types.

isaacsas commented 2 years ago

Could you upload the .net file it generates? That's what we'd want to be able to load here.

pascalschulthess commented 2 years ago

There you go... just remove the .txt suffix. Gerosa_BioNetGen_model.net.txt

isaacsas commented 2 years ago

Thanks! This looks like we could probably support it in a fairly straightforward manner since the observables are never used in reactions or rate laws. The latter would be more complicated to support. The groups should already work (i.e. you should be able to index a solution object using the group names to get their values).

pascalschulthess commented 2 years ago

Alright. I just discovered that if I do:

model_bngl = loadrxnetwork(BNGNetwork(), model_path)
rn = model_bngl.rn

I can access the observables if I do rn.observed which gives for instance tDUSP(t) ~ S24(t) + S296(t). I, however, don't know how I could get this as an output of solve.

isaacsas commented 2 years ago

Here is an example on how to index by the name of a variable (which should work for observables too)

using Catalyst, OrdinaryDiffEq
rn = @reaction_network begin
    α, S + I --> 2I
    β, I --> R
end α β
p     = [:α => .1/1000, :β => .01]
tspan = (0.0,250.0)
u0    = [:S => 999.0, :I => 1.0, :R => 0.0]
op    = ODEProblem(rn, u0, tspan, p)
sol   = solve(op, Tsit5()) 
@unpack S = rn
sol[S]   # solution for the variable S
# or
t = [1.0, 1.5, 2.0]
@unpack I = rn
sol(t, idxs=[S,I])   # values of S and I at the times