Constraints, flattening, and extending for ReactionSystems

[isaacsas]

Enables extending ReactionSystems with algebraic and ODE constraints. While we currently can have sub-systems of algebraic (or even ODE) constraints, this is annoying for specifying algebraic equations to couple systems (since everything needs to get parent scoped). i.e. imagine wanting to hook three systems representing genes together, conceptually this is not necessarily a hierarchical system, but more a collection of parallel systems with coupling equations one level up. (This wouldn't necessarily be an issue if we could make the top-level system a NonlinearSystem, but currently MT requires sub-systems to have the same type, so top-level systems must be ReactionSystems.)

This PR allows attaching to a ReactionSystem an internal constraint system, which is considered at the same scope as the ReactionSystem. In this way we can logically extend ReactionSystems with equations and/or other ReactionSystems (we already support compose). For a simple alias_elimination of basic A ~ B type equations we could then flatten the full hierarchy of systems into one ReactionSystem and one internal sys.constraint system, and then apply the elimination.

Still TODO:

• get_constraints(network) and support a parallel constraint subsystem
• flatten
• extend
• tests
• Allow / use ODESystems as constraints, and error if then try to convert to NonlinearSystem.
• Allow interpolation of variables storing network names in the DSL.

[isaacsas]

Some simple examples

rn = @reaction_network rn begin
    (k1,k2), A <--> B
    end k1 k2
@parameters a,b
@variables t,A(t),C(t)
D = Differential(t)
eqs = [D(C) ~ -b*C + a*A]
@named osys = ODESystem(eqs,t,[A,C],[a,b])
rn2 = extend(osys,rn)
rnodes = convert(ODESystem,rn2)

# if there are actually ODEs in the ODESystem we error
@test_throws ErrorException convert(NonlinearSystem,rn2)

# no ODEs below, so can represent as ODE or Nonlinear systems and convert to either
eqs = [0 ~ -a*A + C, 0 ~ -b*C + a*A]
@named nlsys = NonlinearSystem(eqs,[A,C],[a,b])
rn2 = extend(nlsys,rn)
rnodes = convert(ODESystem,rn2)
rnnlsys = convert(NonlinearSystem,rn2)

@named nlsys = ODESystem(eqs,t,[A,C],[a,b])
rn2 = extend(nlsys,rn)
rnodes = convert(ODESystem,rn2)
rnnlsys = convert(NonlinearSystem,rn2)

[isaacsas]

With flatten we should now have most of the setup to tackle a version of alias_elimination on a ReactionSystem. The goal here is to be able to completely eliminate algebraic coupling equations and spit out a ReactionSystem with no constraints (but some observables). This would really facilitate compositional modeling of networks since such a pure ReactionSystem can then be used in standard network analysis tooling.

I'll probably tackle the alias_elimination in two steps, first for simple A ~ B type equations, and then more complicated linear equations. Supporting the latter would allow us to reduce out conservation laws I think, but is much more complicated. It'd be nice to reuse some of the existing code in MT if I can figure it out easily...

[isaacsas]

Perhaps a better example, which I plan to put in a tutorial:

function repressed_gene(; name)
  @parameters t α₀ α K n δ β μ
  @variables m(t) P(t) R(t)
  rxs = [
         Reaction(α₀, nothing, [m]),
         Reaction(α / (1 + (R/K)^n), nothing, [m]),
         Reaction(δ, [m], nothing),
         Reaction(β, [m], [m,P]),
         Reaction(μ, [P], nothing)
        ]
  
  specs = [m,P,R]
  pars  = [α₀,α,K,n,δ,β,μ]
  ReactionSystem(rxs, t, specs, pars; name=name)
end
@named sys₁ = repressed_gene()
@named sys₂ = repressed_gene()
@named sys₃ = repressed_gene()
connections = [sys₁.R ~ sys₃.P,
               sys₂.R ~ sys₁.P,
               sys₃.R ~ sys₂.P]
@named csys = NonlinearSystem(connections, [sys₁.R,sys₃.P,sys₂.R,sys₁.P,sys₃.R,sys₂.P],[])
@named repressilator = ReactionSystem(t; constraints=csys, systems=[sys₁,sys₂,sys₃])

# or 
@named repressilator = ReactionSystem(t; systems=[sys₁,sys₂,sys₃])
@named repressilator = extend(csys, repressilator)

[TorkelE]

The repressilator example is really elegant! Have been lagging a bit on developments in this direction, but looking forward to playing around with it.

[isaacsas]

The repressilator example is really elegant! Have been lagging a bit on developments in this direction, but looking forward to playing around with it.

Thanks! BTW, do you know if there is a way we could make the macro work in such a function? i.e. something like

function repressed_gene(; name)
    @reaction_network $name begin
    ...
    end
end

It would be cool if we could use the DSL to generate the core gene module in this way.

[TorkelE]

I think it should be possible, but uncertain how to get it working. E.g. I tried

name = :repressed_gene
rg =  @reaction_network :($(name)) begin
    (p,d), 0 <--> X
end p d
rg.name

but got an error

LoadError: AssertionError: @parameters expects a tuple of expressions or an expression of a tuple (`@parameters x y z(t) v[1:3] w[1:2,1:4]` or `@parameters x y z(t) v[1:3] w[1:2,1:4] k=1.0`)
in expression starting at /home/torkelloman/.julia/packages/Catalyst/RTdmC/src/reaction_network.jl:178

Stacktrace:
 [1] _parse_vars(macroname::Symbol, type::Type, x::Tuple{Symbol, Expr, Symbol, Symbol}, transform::Function)
   @ Symbolics ~/.julia/packages/Symbolics/sITWZ/src/variable.jl:134
 [2] var"@parameters"(__source__::LineNumberNode, __module__::Module, xs::Vararg{Any, N} where N)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/Mo4gw/src/parameters.jl:30
 [3] eval
   @ ./boot.jl:360 [inlined]
 [4] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1094

there are some meta-programming stuff going on here which I don't fully remember how to deal with. It would probably be possible to make a function remake(::ReactionSystem;kwargs...) where one could set remake(repressed_gene_base,name=my_name), but that doesn't really solve the initial problem.

[isaacsas]

@TorkelE I think I just added that feature in. Could you take a look at the changes in reaction_network.jl? This seems to handle it and I think I haven't broken anything after trying the various permutations of with/without names that we allow.

[TorkelE]

I think that looks good. Admittedly I haven't done metaprogramming in a while, so a bit rusty when it comes to these things. I'll play around with things a bit when it is merged as well in case I were to find anything, but I'm not worried.

[isaacsas]

Once this gets merged I'm planning to add a tutorial using all this stuff.