@add_reactions macro issue

[isaacsas]

See https://discourse.julialang.org/t/proper-syntax-for-adding-reactions-to-a-reaction-network-using-addreaction-function-in-catalyst-jl/44804

Trying that example (with the parameters added at the end), when doing:

@add_reactions coupled_switch begin
       gD*(Hsn(tBD,B,fBD,nBD)*Hsn(tCD,C,fCD,nCD)), ∅ → D 
       end gD tBD fBD nBD tCD fCD nCD

I get an error that B is not defined. @TorkelE does the macro not allow using species that were previously defined in the network within rates?

[ChrisRackauckas]

I think it might be colliding a bit with https://github.com/SciML/Catalyst.jl/issues/248 as well. I think there's something funny in the way the functions are handled that's causing this to fail, but I'm not quite sure.

[isaacsas]

I think it is separate as this works:

@add_reactions coupled_switch begin
       gD, ∅ → D 
end gD

but this fails

@add_reactions coupled_switch begin
       gD*B, ∅ → D 
end gD

[isaacsas]

This also works

@add_reactions coupled_switch begin
       gD*D, ∅ → D 
end gD

I think we need to update the macro to check for species that are in coupled_switch as allowed symbols.

[ChrisRackauckas]

ahh I see

[TorkelE]

Yes, that's it. Since what the @add_reactions macro essentially does is creating a new ReactionSystem of your network, and then merge it with the previous one, every set of input reactions that would give an error if following @reaction_network will give the same error if following @add_reactions.

I am not sure why, but I think we decided that every parameter declared in a reaction network have to be listed at the end, and every species must be mentioned in a reaction. Otherwise, there's an error, which is what happens here. This probably makes sense. Either we change this, or we have to implement a custom macro for @add_reactions which can take this into account.

[isaacsas]

Hmm, does it make sense to make anything that doesn't appear in the parameter list or in a reaction a species by default? Maybe that could also useful in building coupled systems? Merging should make sure to not double include a species, so it shouldn't be problematic to have this.

[TorkelE]

I vaguely think this was the case in the beginning, but it is not that any longer.

The disadvantage of the current system is that it is less flexible. Also, changing it would allow this by default. I think the advantage is that it might help to catch some otherwise hard to find bugs, e.g. if you declare a model

@reaction_network begin
       k1*X1, ∅ → X1
       k2*X2, X1 → X2
       k3*X3, X2 → X3
       k4*x4, X3 → X4
       k5*X5, X4 → X5
       k6*X6, X6 → ∅
end k1 k2 k3 k4 k5 k6

but here mistyped a X4 as x4, then this currently this yields an error (since species most likely exists in a reaction). Updating it will let this model parse, and it might take a while until the bug is discovered. (If we were to rewrite the @add_reactions macro, we could make a new @reaction_network macro, which does parse these networks, and can be used in cases when the user know theres an exception)

[isaacsas]

The behavior should probably be consistent between @add_reactions and @reaction_network, so either we always allow undefined symbols to become species, we add a parallel set of macros that assumes this, or we have another way to handle it. Could we add a way to designate something is a species within the parameter list? i.e. like

@reaction_network begin
       k1*X1, ∅ → X1
       k2*X2, X1 → X2
       k3*X3, X2 → X3
       k4*A*B, X3 → X4
       k5*X5, X4 → X5
       k6*X6, X6 → ∅
end k1 k2 k3 k4 k5 k6, A B

Here A and B are added as species. (I have no idea what, if any notation is reasonable and will work in a macro here.) I do like the current behavior of giving an error since it is useful for finding bugs.

[isaacsas]

Otherwise I guess we just modify the docs to tell people to add fictitious reactions?

@reaction_network begin
       k1*X1, ∅ → X1
       k2*X2, X1 → X2
       k3*X3, X2 → X3
       k4*A*B, X3 → X4
       k5*X5, X4 → X5
       k6*X6, X6 → ∅
       0.0, A --> A
       0.0, B --> B
end k1 k2 k3 k4 k5 k6, A B

Seems a bit hacky, but it should work I think (I haven't tested it beyond a simple SSA example though).

[TorkelE]

The second one seems good for now.

If this is something that will actually happen occasionally, we probably should implement a less hacky solution.

@reaction_network begin
       k1*X1, ∅ → X1
       k2*X2, X1 → X2
       k3*X3, X2 → X3
       k4*A*B, X3 → X4
       k5*X5, X4 → X5
       k6*X6, X6 → ∅
end k1 k2 k3 k4 k5 k6, A B

seems like a good idea.

For the @add_reactions, I am personally fine with letting it work a little bit different. E.g. I would be fine with not requiring a parameter to be declared if it was already declared in the base network. Same allowing species without declaring them as such, as long as they were somehow declared in the initial network.

[isaacsas]

Yeah, that is reasonable too. I guess we don't need to require the same behavior since it is an extension macro. If something is in the base network allow it to be used in @add_reactions too.

[TorkelE]

I'd say that if

rn = @reaction_network begin
    r1s...
end p1...

@add_reactions rn begin
    r2s...
end p2...

is a valid notation, then

rn = @reaction_network begin
    r1s...
    r2s...
end p1... p2...

should be, but not necessarily

rn = @reaction_network begin
    r2s...
end p2...

[isaacsas]

This original issue seems resolved on master using @add_reactions so going to close.