Species' order is not defined by @species

[johannesnauta]

I encountered a problem where the order of species is not defined by @species. This lead to problem when creating/defining the ODEProblem, as I expected the order of the initial state to follow the order as in @species. More specifically, consider the following example with 2 species:

julia> using Catalyst

julia> @parameters α[1:2]
1-element Vector{Symbolics.Arr{Num, 1}}:
 α[1:2]

julia> @variables t
1-element Vector{Num}:
 t

julia> @species (x(t))[1:2]
1-element Vector{Symbolics.Arr{Num, 1}}:
 (x(t))[1:2]

julia> r1 = Reaction(α[1], [x[1]], [x[1]], [1], [2])
α[1], (x(t))[1] --> 2*(x(t))[1]

julia> r2 = Reaction(α[2], [x[2]], [x[2]], [1], [2])
α[2], (x(t))[2] --> 2*(x(t))[2]

julia> @named rs1 = ReactionSystem([r1,r2])
Model rs1
States (2):
  (x(t))[1]
  (x(t))[2]
Parameters (2):
  α[1]
  α[2]

julia> @named rs2 = ReactionSystem([r2,r1])
Model rs2
States (2):
  (x(t))[2]
  (x(t))[1]
Parameters (2):
  α[2]
  α[1]

Note the symbolic definition as the systems I am working with are in general of some length n. In the example, I create the growth/replication reactions with some rate α. When creating the ODEProblem from the reaction you need to provide initial conditions u0, but these of course depend on the actual order of species(rs). As I expected the species to be ordered as [x(t)[1], x(t)[2]] the solutions to the system of ODEs was nonsensical.

Is this behavior intended? Is there a way to pre-define the order of the species using @species x[1:n] (or something similar)?

[johannesnauta]

I have noticed that one can specify the order of species when defining the ReactionSystem, as shown in the example in the documentation. However, in order to do this, I need to also provide the parameters of the system (called ps in the source code).
The issue now is that I do not have a single parameter, as k[1:n] in the example. What I have is more like:

@parameters α[1:n,1:n], η[1:n]
@variables t
@species (x(t))[1:n]

Note the fact that α is a matrix, such as e.g. in the Lotka-Volterra system where α encodes the competition. I have tried multiple ways to provide ps, that all fail:

ps = [α, η]  #Raises ERROR: MethodError: no method matching nameof(::Vector{Symbolics.Num})
ps = vcat([collect(z) for z in [α, η]])  #Raises ERROR: MethodError: no method matching nameof(::Vector{Symbolics.Num})

Is there a way to provide the (order of) species without providing the parameters? Or is there a way multiple parameters should be presented in order for the ReactionSystem to be created without failing?

[TorkelE]

@isaacsas Is the one of us most familiar with creating ReactionSystems directly using the components, while I usually use the macro. I think there should be a way to reorder them while using your approach,

The @species command only sets the order of the species when used within the macro. E.g.

using Catalyst
rs = @reaction_network begin
    @species Y(t) X(t)
    p1, ∅ → X
    p2, ∅ → Y
end
species(rs)

returns [Y(t), X(t)]. When it is used outside of the macro, @species simply designates that these symbols denote species, but do not influence their order within the ReactionSystems within which they are used.

[johannesnauta]

Makes sense. I was just confused as I did not anticipate the order change before. For now, I just change the order of the reactions when defining the ReactionSystem, but this is not robust. Looking through the source code it appears you can define the species order, but you have to provide the parameters in this case as well. I am not sure why there is no option to give one of them, but there is most likely some underlying reason why this is not straightforward to do as I expect that the parameters might need some specific order as well in this case.

[TorkelE]

I agree that setting these independently would be useful.

Also, I think it is a general policy of https://github.com/SciML/ModelingToolkit.jl that one should not depend on the order of variables and parameters (which then transfers to Catalyst). Generally, you should be able to use e.g X as an index instead of 1. Is this not possible in your case?

[johannesnauta]

Generally, you should be able to use e.g X as an index instead of 1. Is this not possible in your case?

I am not sure what you mean by this. If I have n species, then I need to define them as @species x[1:n], right? Then you can index them with i of course.

[TorkelE]

E.g here you can use the species to fetch what you want from e.g. a solution object:

@variables t
@species x(t)[1:2] y(t)[1:2]
r1 = Reaction(1, [x[1]], [x[2]], [1], [2])
r2 = Reaction(1, [y[1]], [y[2]], [1], [3])
@named rs1 = ReactionSystem([r1, r2])
oprob = ODEProblem(rs1, [1.0, 0.0, 2.0, 0.0], (0.0,10.0))
sol = solve(oprob, Tsit5())

sol[x[2]]
sol[y[1]]

if that make sense?

[johannesnauta]

Ah that makes sense. I did not know that you could index the solution with the species -- now I see what you mean. This is indeed another good way of doing things. I do not need to preserve the order as long as I use the index in the solution. Thanks for that insight!

[isaacsas]

This will give you a vector of all the parameters you should be able to use with ReactionSystem:

ps = vcat([vec(collect(z)) for z in [α, η]]...)

vcat requires each vector to concatenate to be a separate argument (i.e. it doesn't concatenate a vector of vectors as a single argument), i.e. vcat(v1, v2, v3) not vcat([v1,v2,v3]).

or just loop and append! them into a new vector:

ps_to_merge = @parameters α[1:n,1:n], η[1:n], a
ps = Num[]
for p in ps_to_merge
    pv = vec(collect(p))
    append!(ps, pv)
end

[TorkelE]

This gives some examples of places where you might want to inference with stuff like problems: https://docs.sciml.ai/Catalyst/dev/catalyst_applications/simulation_structure_interfacing/

If you find any edge cases where you cannot use species, do tell us!

[isaacsas]

@johannesnauta I'm going to close as I think we've answered your questions, but feel free to reopen if there are further problems.

[johannesnauta]

This gives some examples of places where you might want to inference with stuff like problems: https://docs.sciml.ai/Catalyst/dev/catalyst_applications/simulation_structure_interfacing/

I knew there must have been a way to define the initial state u0 using (Symbolic) pairs, but somehow I could not find the documentation that you just provided. I will keep that page in mind when I encounter similar things in the future.

[isaacsas]

I completely forgot probably the easiest way:

reduce(vcat, vec(collect(z)) for z in [α, η])

I should have said, the call to vec is needed as collect returns a matrix and not a vector when applied to a matrix.

[isaacsas]

See also the first tutorial's section on the repressilator, which shows how to use symbolic pairs for the parameter and initial condition maps:

https://docs.sciml.ai/Catalyst/dev/introduction_to_catalyst/introduction_to_catalyst/#Mass-action-ODE-models

[johannesnauta]

Yeah I should've used the u0map way before, this would probably not have lead me to prompting this question even. I just falsely assumed what the @species would do without looking into it. Took me quite some time to see that it just switched around the order of the state variables.

[isaacsas]

You can also specify an initial condition / parameter value vector/matrix using default values in @species and @parameters. i.e.

using Catalyst
n = 2
@parameters α[1:n,1:n]=rand(n,n) η[1:n]=rand(n)
@variables t
@species (x(t))[1:n]=rand(n)
@named rx = ReactionSystem(Reaction[], t, collect(x), vcat(vec(collect(α)), collect(η)))

giving

Model rx
States (2):
  (x(t))[1] [defaults to 0.175127]
  (x(t))[2] [defaults to 0.477996]
Parameters (6):
  α[1, 1] [defaults to 0.616545]
  α[2, 1] [defaults to 0.121818]
  α[1, 2] [defaults to 0.625685]
  α[2, 2] [defaults to 0.93927]
  η[1] [defaults to 0.437466]
  η[2] [defaults to 0.0977467]

[johannesnauta]

That might also be useful, at least for the species (in my case). Thanks for all the help!