Documentation doesn't explain how to programmatically generate species

[nathanielvirgo]

The documentation at https://docs.sciml.ai/Catalyst/stable/catalyst_functionality/programmatic_CRN_construction/ explains how to generate reactions programmatically, but it doesn't mention how to do this for species. It would obviously be a huge drawback if you had to enter the list of species manually each time!

I'm submitting this as a bug report because I think this should be fixed in the documentation, but in addition I would also just like to know how to do this.

for example, what should I do if I want to generate this reaction network?

A1 <=> A2 A2 <=> A3 ... A999 <=> A1000

[isaacsas]

Does this example help?

https://docs.sciml.ai/Catalyst/stable/example_networks/smoluchowski_coagulation_equation/

[nathanielvirgo]

It does help, but in terms of documentation it's not very clear - it's a very complicated example and that particular aspect of it isn't really explained. I guess the relevant line is

@species k[1:nr] (X(t))[1:N]

and I guess the k[1:nr] part is optional, since that seems to be declaring reaction rates and not species. (Which seems kind of weird - why aren't they declared with @parameters instead?) I can extrapolate from this by guesswork but it would be nice to have it documented properly.

It would also be helpful to know how to do the following things, if it's possible to do them:

• add a new species to the network dynamically - for some algorithms it might make sense to add a few reactions, then add a few species, then more reactions and so on, rather than declaring all the species at the start
• give the species names other than X[i]. In the example I gave calling them X[i] makes sense, but it's easy to imagine cases where you'd want to name them differently, with the names also generated programmatically in the form of strings.

[isaacsas]

Yeah, k should be a parameter. It is a relic from a much older version of Catalyst where there was less of a distinction, and looks like we missed updating it.

You can create a species from a string as follows:

str = "AB"
strsym = Symbol(str)
@variables t
spcs = @species $(strsym)(t)

giving

julia> spcs[1]
AB(t)

We should indeed update the programmatic tutorial to show this.

If you want to add new species and/or reactions I'd suggest creating a new network using extend and/or compose as shown at https://docs.sciml.ai/Catalyst/stable/catalyst_functionality/compositional_modeling/

These take one or more systems and merge them together / create a tree-like system from them.

Alternatively, we do provide addspecies!, addparam! and addreaction! to mutate a system, see

https://docs.sciml.ai/Catalyst/stable/api/catalyst_api/#Functions-to-extend-or-modify-a-network

however, I'd recommend not using these unless really needed as they may be deprecated in the near future. (The reason being that ModelingToolkit wants systems to be immutable, with adding/removing components done via generating new systems. So we have debated removing these features to maintain compatibility with that philosophy, and don't really advertise them in the docs.)

[nathanielvirgo]

Great, that's really helpful, thank you

[nathanielvirgo]

I'm having trouble using species created this way in reactions:

@variables t
@parameters α

# species created normally
@species x(t)
@species y(t)

# species created programmatically
strsym = Symbol("A")
a = @species $(strsym)(t)
strsym = Symbol("B")
b = @species $(strsym)(t)

# 'normal' species work ok:
Reaction(α,[x],[y])

# but programmatically created ones do not:
Reaction(α,[a],[b])

gives the error

MethodError: no method matching getmetadata(::Vector{Num}, ::Type{Catalyst.VariableSpecies}, ::Bool)

Closest candidates are:
  getmetadata(::SymbolicUtils.Symbolic, ::Any, ::Any)
   @ SymbolicUtils ~/.julia/packages/SymbolicUtils/Oyu8Z/src/types.jl:629
  getmetadata(::SymbolicUtils.Symbolic, ::Any)
   @ SymbolicUtils ~/.julia/packages/SymbolicUtils/Oyu8Z/src/types.jl:620
  getmetadata(::Complex{Num}, ::Any)
   @ Symbolics ~/.julia/packages/Symbolics/BQlmn/src/Symbolics.jl:156
  ...


Stacktrace:
  [1] isspecies(s::Vector{Num})
    @ Catalyst ~/.julia/packages/Catalyst/zTvWJ/src/reactionsystem.jl:36
  [2] isvalidreactant(s::Vector{Num})
    @ Catalyst ~/.julia/packages/Catalyst/zTvWJ/src/reactionsystem.jl:126
  [3] _all
    @ ./reduce.jl:1283 [inlined]
  [4] #all#830
    @ ./reducedim.jl:1007 [inlined]
  [5] all
    @ ./reducedim.jl:1007 [inlined]
  [6] Reaction(rate::Num, subs::Vector{Vector{Num}}, prods::Vector{Vector{Num}}, substoich::Vector{Int64}, prodstoich::Vector{Int64}; netstoich::Nothing, only_use_rate::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Catalyst ~/.julia/packages/Catalyst/zTvWJ/src/reactionsystem.jl:172
  [7] Reaction
    @ ~/.julia/packages/Catalyst/zTvWJ/src/reactionsystem.jl:128 [inlined]
  [8] Reaction(rate::Num, subs::Vector{Vector{Num}}, prods::Vector{Vector{Num}}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Catalyst ~/.julia/packages/Catalyst/zTvWJ/src/reactionsystem.jl:191
  [9] Reaction(rate::Num, subs::Vector{Vector{Num}}, prods::Vector{Vector{Num}})
    @ Catalyst ~/.julia/packages/Catalyst/zTvWJ/src/reactionsystem.jl:188
 [10] top-level scope
    @ In[63]:18

[isaacsas]

See my example above; @species returns a vector with the first entry being your symbolic variable.