The Reaction DSL
This tutorial describes the syntax for building chemical reaction network models using Catalyst's domain-specific language (DSL). Examples showing how to both construct and solve ODE, SDE, and jump models are provided in Basic Chemical Reaction Network Examples. To learn more about the symbolic ReactionSystem
s generated by the DSL, and how to use them directly, see the tutorial on Programmatic Construction of Symbolic Reaction Systems.
We first load the Catalyst
package, which is required for the code in this tutorial to run
using Catalyst
Basic syntax
The @reaction_network
macro allows the (symbolic) specification of reaction networks with a simple format. Its input is a set of chemical reactions, and from them it generates a symbolic ReactionSystem
reaction network object. The ReactionSystem
can be used as input to ModelingToolkit ODEProblem
, NonlinearProblem
, SteadyStateProblem
, SDEProblem
, JumpProblem
, and more. ReactionSystem
s can also be incrementally extended as needed, allowing for programmatic construction of networks and network composition.
The basic syntax is:
rn = @reaction_network begin
2.0, X + Y --> XY
1.0, XY --> Z1 + Z2
end
\[ \begin{align*} \mathrm{X} + \mathrm{Y} &\xrightarrow{2.0} \mathrm{XY} \\ \mathrm{XY} &\xrightarrow{1.0} \mathrm{Z1} + \mathrm{Z2} \end{align*} \]
where each line of the @reaction_network
macro corresponds to a chemical reaction. Each reaction consists of a reaction rate (the expression on the left-hand side of ,
), a set of substrates (the expression in-between ,
and -->
), and a set of products (the expression on the right-hand side of -->
). The substrates and the products may contain one or more reactants, separated by +
. The naming convention for these is the same as for normal variables in Julia.
The chemical reaction model is generated by the @reaction_network
macro and stored in the rn
variable (a normal Julia variable, which does not need to be called rn
). It corresponds to a ReactionSystem
, a symbolic representation of the chemical network. The generated ReactionSystem
can be converted to a symbolic differential equation model via
osys = convert(ODESystem, rn)
\[ \begin{align} \frac{\mathrm{d} X\left( t \right)}{\mathrm{d}t} =& - 2 Y\left( t \right) X\left( t \right) \\ \frac{\mathrm{d} Y\left( t \right)}{\mathrm{d}t} =& - 2 Y\left( t \right) X\left( t \right) \\ \frac{\mathrm{d} \mathrm{XY}\left( t \right)}{\mathrm{d}t} =& - \mathrm{XY}\left( t \right) + 2 Y\left( t \right) X\left( t \right) \\ \frac{\mathrm{d} \mathrm{Z1}\left( t \right)}{\mathrm{d}t} =& \mathrm{XY}\left( t \right) \\ \frac{\mathrm{d} \mathrm{Z2}\left( t \right)}{\mathrm{d}t} =& \mathrm{XY}\left( t \right) \end{align} \]
We can then convert the symbolic ODE model into a compiled, optimized representation for use in the SciML ODE solvers by constructing an ODEProblem
. Creating an ODEProblem
also requires our specifying the initial conditions for the model. We do this by creating a mapping from each symbolic variable representing a chemical species to its initial value
# define the symbolic variables
@variables t
@species X(t) Y(t) Z(t) XY(t) Z1(t) Z2(t)
# create the mapping
u0 = [X => 1.0, Y => 1.0, XY => 1.0, Z1 => 1.0, Z2 => 1.0]
5-element Vector{Pair{Num, Float64}}:
X(t) => 1.0
Y(t) => 1.0
XY(t) => 1.0
Z1(t) => 1.0
Z2(t) => 1.0
Alternatively, we can create a mapping using Julia Symbol
s for each variable, and then convert them to a mapping involving symbolic variables like
u0 = symmap_to_varmap(rn, [:X => 1.0, :Y => 1.0, :XY => 1.0, :Z1 => 1.0, :Z2 => 1.0])
5-element Vector{Pair{Num, Float64}}:
X(t) => 1.0
Y(t) => 1.0
XY(t) => 1.0
Z1(t) => 1.0
Z2(t) => 1.0
Given the mapping, we can then create an ODEProblem
from our symbolic ODESystem
tspan = (0.0, 1.0) # the time interval to solve on
oprob = ODEProblem(osys, u0, tspan, [])
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1.0)
u0: 5-element Vector{Float64}:
1.0
1.0
1.0
1.0
1.0
Catalyst provides a shortcut to avoid having to explicitly convert
to an ODESystem
and/or use symmap_to_varmap
, allowing direct construction of the ODEProblem
like
u0 = [:X => 1.0, :Y => 1.0, :XY => 1.0, :Z1 => 1.0, :Z2 => 1.0]
oprob = ODEProblem(rn, u0, tspan, [])
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1.0)
u0: 5-element Vector{Float64}:
1.0
1.0
1.0
1.0
1.0
For more detailed examples, see the Basic Chemical Reaction Network Examples.
Defining parameters and species
Numeric parameter values do not need to be set when the model is created, i.e. Catalyst supports symbolic parameters too:
rn = @reaction_network begin
k1, X --> Y
k2, Y --> X
end
\[ \begin{align*} \mathrm{X} &\xrightleftharpoons[k2]{k1} \mathrm{Y} \end{align*} \]
All symbols that do not appear as a substrate or product in a reaction are designated by Catalyst as a parameter (i.e. all symbols appearing only within rate expressions and/or as stoichiometric coefficients). In this example X
and Y
appear as a substrates and products, but neither k1
nor k2
. Hence k1
and k2
are designated as parameters. Later in this tutorial, we will describe how to manually specify what should be considered a species or parameter.
Production, Destruction, and Stoichiometry
Sometimes reactants are produced/destroyed from/to nothing. This can be designated using either 0
or ∅
:
rn = @reaction_network begin
2.0, 0 --> X
1.0, X --> 0
end
\[ \begin{align*} \varnothing &\xrightleftharpoons[1.0]{2.0} \mathrm{X} \end{align*} \]
If several molecules of the same reactant are involved in a reaction, the stoichiometry of a reactant in a reaction can be set using a number. Here, two molecules of species X
form the dimer X2
:
rn = @reaction_network begin
1.0, 2X --> Y
end
\[ \begin{align*} 2 \mathrm{X} &\xrightarrow{1.0} \mathrm{Y} \end{align*} \]
this corresponds to the differential equation:
convert(ODESystem, rn)
\[ \begin{align} \frac{\mathrm{d} X\left( t \right)}{\mathrm{d}t} =& - \left( X\left( t \right) \right)^{2} \\ \frac{\mathrm{d} Y\left( t \right)}{\mathrm{d}t} =& \frac{1}{2} \left( X\left( t \right) \right)^{2} \end{align} \]
Other numbers than 2 can be used, and parenthesis can be used to reuse the same stoichiometry for several reactants:
rn = @reaction_network begin
1.0, X + 2(Y + Z) --> W
end
\[ \begin{align*} \mathrm{X} + 2 \mathrm{Y} + 2 \mathrm{Z} &\xrightarrow{1.0} \mathrm{W} \end{align*} \]
Note, one can explicitly multiply by integer coefficients too, i.e.
rn = @reaction_network begin
1.0, X + 2*(Y + Z) --> W
end
\[ \begin{align*} \mathrm{X} + 2 \mathrm{Y} + 2 \mathrm{Z} &\xrightarrow{1.0} \mathrm{W} \end{align*} \]
Arrow variants
A variety of Unicode arrows are accepted by the DSL in addition to -->
. All of these work: >
, →
↣
, ↦
, ⇾
, ⟶
, ⟼
, ⥟
, ⥟
, ⇀
, ⇁
. Backwards arrows can also be used to write the reaction in the opposite direction. For example, these reactions are equivalent:
rn = @reaction_network begin
1.0, X + Y --> XY
1.0, X + Y → XY
1.0, XY ← X + Y
1.0, XY <-- X + Y
end
\[ \begin{align*} \mathrm{X} + \mathrm{Y} &\xrightarrow{1.0} \mathrm{XY} \\ \mathrm{X} + \mathrm{Y} &\xrightarrow{1.0} \mathrm{XY} \\ \mathrm{X} + \mathrm{Y} &\xrightarrow{1.0} \mathrm{XY} \\ \mathrm{X} + \mathrm{Y} &\xrightarrow{1.0} \mathrm{XY} \end{align*} \]
Bi-directional arrows for reversible reactions
Bi-directional arrows, including bidirectional Unicode arrows like ↔, can be used to designate a reversible reaction. For example, these two models are equivalent:
rn = @reaction_network begin
2.0, X + Y --> XY
2.0, X + Y <-- XY
end
\[ \begin{align*} \mathrm{X} + \mathrm{Y} &\xrightleftharpoons[2.0]{2.0} \mathrm{XY} \end{align*} \]
rn2 = @reaction_network begin
(2.0,2.0), X + Y <--> XY
end
\[ \begin{align*} \mathrm{X} + \mathrm{Y} &\xrightleftharpoons[2.0]{2.0} \mathrm{XY} \end{align*} \]
If the reaction rates in the backward and forward directions are different, they can be designated in the following way:
rn = @reaction_network begin
(2.0,1.0), X + Y <--> XY
end
\[ \begin{align*} \mathrm{X} + \mathrm{Y} &\xrightleftharpoons[1.0]{2.0} \mathrm{XY} \end{align*} \]
which is identical to
rn = @reaction_network begin
2.0, X + Y --> XY
1.0, X + Y <-- XY
end
\[ \begin{align*} \mathrm{X} + \mathrm{Y} &\xrightleftharpoons[1.0]{2.0} \mathrm{XY} \end{align*} \]
Combining several reactions in one line
Several similar reactions can be combined in one line by providing a tuple of reaction rates and/or substrates and/or products. If several tuples are provided, they must all be of identical length. These pairs of reaction networks are all identical.
Pair 1:
rn1 = @reaction_network begin
1.0, S --> (P1,P2)
end
\[ \begin{align*} \mathrm{S} &\xrightarrow{1.0} \mathrm{P1} \\ \mathrm{S} &\xrightarrow{1.0} \mathrm{P2} \end{align*} \]
rn2 = @reaction_network begin
1.0, S --> P1
1.0, S --> P2
end
\[ \begin{align*} \mathrm{S} &\xrightarrow{1.0} \mathrm{P1} \\ \mathrm{S} &\xrightarrow{1.0} \mathrm{P2} \end{align*} \]
Pair 2:
rn1 = @reaction_network begin
(1.0,2.0), (S1,S2) --> P
end
\[ \begin{align*} \mathrm{S1} &\xrightarrow{1.0} \mathrm{P} \\ \mathrm{S2} &\xrightarrow{2.0} \mathrm{P} \end{align*} \]
rn2 = @reaction_network begin
1.0, S1 --> P
2.0, S2 --> P
end
\[ \begin{align*} \mathrm{S1} &\xrightarrow{1.0} \mathrm{P} \\ \mathrm{S2} &\xrightarrow{2.0} \mathrm{P} \end{align*} \]
Pair 3:
rn1 = @reaction_network begin
(1.0,2.0,3.0), (S1,S2,S3) → (P1,P2,P3)
end
\[ \begin{align*} \mathrm{S1} &\xrightarrow{1.0} \mathrm{P1} \\ \mathrm{S2} &\xrightarrow{2.0} \mathrm{P2} \\ \mathrm{S3} &\xrightarrow{3.0} \mathrm{P3} \end{align*} \]
rn2 = @reaction_network begin
1.0, S1 --> P1
2.0, S2 --> P2
3.0, S3 --> P3
end
\[ \begin{align*} \mathrm{S1} &\xrightarrow{1.0} \mathrm{P1} \\ \mathrm{S2} &\xrightarrow{2.0} \mathrm{P2} \\ \mathrm{S3} &\xrightarrow{3.0} \mathrm{P3} \end{align*} \]
This can also be combined with bi-directional arrows, in which case separate tuples can be provided for the backward and forward reaction rates. These reaction networks are identical
rn1 = @reaction_network begin
(1.0,(1.0,2.0)), S <--> (P1,P2)
end
\[ \begin{align*} \mathrm{S} &\xrightarrow{1.0} \mathrm{P1} \\ \mathrm{S} &\xrightarrow{1.0} \mathrm{P2} \\ \mathrm{P1} &\xrightarrow{1.0} \mathrm{S} \\ \mathrm{P2} &\xrightarrow{2.0} \mathrm{S} \end{align*} \]
rn2 = @reaction_network begin
1.0, S --> P1
1.0, S --> P2
1.0, P1 --> S
2.0, P2 --> S
end
\[ \begin{align*} \mathrm{S} &\xrightarrow{1.0} \mathrm{P1} \\ \mathrm{S} &\xrightarrow{1.0} \mathrm{P2} \\ \mathrm{P1} &\xrightarrow{1.0} \mathrm{S} \\ \mathrm{P2} &\xrightarrow{2.0} \mathrm{S} \end{align*} \]
Variable reaction rates
Reaction rates do not need to be a single parameter or a number, but can also be expressions depending on time or the current amounts of system species (when, for example, one species can activate the production of another). For instance, this is a valid notation:
rn = @reaction_network begin
1.0, X --> ∅
k*X, Y --> ∅
end
\[ \begin{align*} \mathrm{X} &\xrightarrow{1.0} \varnothing \\ \mathrm{Y} &\xrightarrow{X k} \varnothing \end{align*} \]
corresponding to the ODE model
convert(ODESystem,rn)
\[ \begin{align} \frac{\mathrm{d} X\left( t \right)}{\mathrm{d}t} =& - X\left( t \right) \\ \frac{\mathrm{d} Y\left( t \right)}{\mathrm{d}t} =& - k Y\left( t \right) X\left( t \right) \end{align} \]
With respect to the corresponding mass action ODE model, this is actually equivalent to the reaction system
rn = @reaction_network begin
1.0, X --> ∅
k, X + Y --> X
end
\[ \begin{align*} \mathrm{X} &\xrightarrow{1.0} \varnothing \\ \mathrm{X} + \mathrm{Y} &\xrightarrow{k} \mathrm{X} \end{align*} \]
convert(ODESystem,rn)
\[ \begin{align} \frac{\mathrm{d} X\left( t \right)}{\mathrm{d}t} =& - X\left( t \right) \\ \frac{\mathrm{d} Y\left( t \right)}{\mathrm{d}t} =& - k Y\left( t \right) X\left( t \right) \end{align} \]
While the ODE models corresponding to the preceding two reaction systems are identical, in the latter example the Reaction
stored in rn
will be classified as ismassaction
while in the former it will not, which can impact optimizations used in generating JumpSystem
s. For this reason, it is recommended to use the latter representation when possible.
Most expressions and functions are valid reaction rates, e.g.:
using SpecialFunctions
rn = @reaction_network begin
2.0*X^2, 0 --> X + Y
t*gamma(Y), X --> ∅
pi*X/Y, Y → ∅
end
\[ \begin{align*} \varnothing &\xrightarrow{2.0 X^{2}} \mathrm{X} + \mathrm{Y} \\ \mathrm{X} &\xrightarrow{t \Gamma\left( Y \right)} \varnothing \\ \mathrm{Y} &\xrightarrow{\frac{\pi X}{Y}} \varnothing \end{align*} \]
where here t
always denotes Catalyst's time variable. Please note that many user-defined functions can be called directly, but others will require registration with Symbolics.jl (see the faq).
Explicit specification of network species and parameters
Recall that the @reaction_network
macro automatically designates symbols used in the macro as either parameters or species, with symbols that appear as a substrate or product being species, and all other symbols becoming parameters (i.e. those that only appear within a rate expression and/or as stoichiometric coefficients). Sometimes, one might want to manually override this default behavior for a given symbol. E.g one might want something to be considered as a species, even if it only appears within a rate expression. In the following network
rn = @reaction_network begin
k*X, Y --> 0
end
\[ \begin{align*} \mathrm{Y} &\xrightarrow{X k} \varnothing \end{align*} \]
X
(as well as k
) will be considered a parameter.
By using the @species
and @parameters
options within the @reaction_network
macro, one can manually declare that specified symbols should be considered a species or parameter. E.g in:
rn = @reaction_network begin
@species X(t) Y(t)
k*X, Y --> 0
end
\[ \begin{align*} \mathrm{Y} &\xrightarrow{X k} \varnothing \end{align*} \]
X
and Y
are set as species. Please note that when declaring species using the @species
option, their dependant variable (almost always t
) also needs to be designated. Similarly in
rn = @reaction_network begin
@parameters k
k*X, Y --> 0
end
\[ \begin{align*} \mathrm{Y} &\xrightarrow{X k} \varnothing \end{align*} \]
both X
and k
will be considered as parameters. It is also possible to use both options simultaneously, allowing users to fully specify which symbols are species and/or parameters:
rn = @reaction_network begin
@species X(t) Y(t)
@parameters k
k*X, Y --> 0
end
\[ \begin{align*} \mathrm{Y} &\xrightarrow{X k} \varnothing \end{align*} \]
Here, X
and Y
are designated as species and k
as a parameter.
The lists provided to the @species
and @parameters
options do not need to be extensive. Any symbol that appears in neither list will use the default option as determined by the macro. E.g. in the previous example, where we only want to change the default designation of X
(making it a species rather than a parameter), we can simply write:
rn = @reaction_network begin
@species X(t)
k*X, Y --> 0
end
\[ \begin{align*} \mathrm{Y} &\xrightarrow{X k} \varnothing \end{align*} \]
Finally, note that the @species
and @parameters
options can also be used in begin ... end
block form, allowing more formatted lists of species/parameters:
rn = @reaction_network begin
@parameters begin
d1
d2
end
@species begin
X1(t)
X2(t)
end
d2, X2 --> 0
d1, X1 --> 0
end
\[ \begin{align*} \mathrm{X2} &\xrightarrow{d2} \varnothing \\ \mathrm{X1} &\xrightarrow{d1} \varnothing \end{align*} \]
This can be especially useful when declaring default values for clarity of model specification (see the next section).
Setting default values for initial conditions and parameters
When using the @species
and @parameters
macros to declare species and/or parameters, one can also provide default initial conditions for each species and values for each parameter:
rn = @reaction_network begin
@species X(t)=1.0
@parameters p=1.0 d=0.1
p, 0 --> X
d, X --> ∅
end
\[ \begin{align*} \varnothing &\xrightleftharpoons[d]{p} \mathrm{X} \end{align*} \]
This system can now be simulated without providing initial condition or parameter vectors to the DifferentialEquations.jl solvers:
using DifferentialEquations, Plots
u0 = []
tspan = (0.0, 10.0)
p = []
oprob = ODEProblem(rn, u0, tspan, p)
sol = solve(oprob)
plot(sol)
When providing default values, it is possible to do so for only a subset of the species or parameters, in which case the rest can be specified when constructing the problem type to solve:
rn = @reaction_network begin
@species X(t)
@parameters p=1.0 d
p, 0 --> X
d, X --> 0
end
u0 = [:X => 1.0]
tspan = (0.0, 10.0)
p = [:d => .1]
oprob = ODEProblem(rn, u0, tspan, p)
sol = solve(oprob)
plot(sol)
Finally, default values can be overridden by passing mapping vectors to the DifferentialEquations.jl problem being constructed. Only those initial conditions or parameters for which we want to change their value from the default will need to be passed
u0 = [:X => 1.0]
tspan = (0.0, 10.0)
p = [:p => 2.0, :d => .1] # we change p to 2.0
oprob = ODEProblem(rn, u0, tspan, p)
sol = solve(oprob)
plot(sol)
Setting initial conditions that depend on parameters
It is possible to set the initial condition of one (or several) species so that they depend on some system parameter. This is done in a similar way as default initial conditions, but giving the parameter instead of a value. When doing this, we also need to ensure that the initial condition parameter is a variable of the system:
rn = @reaction_network begin
@parameters X0
@species X(t)=X0
p, 0 --> X
d, X --> ∅
end
\[ \begin{align*} \varnothing &\xrightleftharpoons[d]{p} \mathrm{X} \end{align*} \]
We can now simulate the network without providing any initial conditions:
u0 = []
tspan = (0.0, 10.0)
p = [:p => 2.0, :d => .1, :X0 => 1.0]
oprob = ODEProblem(rn, u0, tspan, p)
sol = solve(oprob)
plot(sol)
Naming the generated ReactionSystem
ModelingToolkit uses system names to allow for compositional and hierarchical models. To specify a name for the generated ReactionSystem
via the @reaction_network
macro, just place the name before begin
:
rn = @reaction_network production_degradation begin
p, ∅ --> X
d, X --> ∅
end
ModelingToolkit.nameof(rn) == :production_degradation
true
Pre-defined functions
Hill functions and a Michaelis-Menten function are pre-defined and can be used as rate laws. Below, the pair of reactions within rn1
are equivalent, as are the pair of reactions within rn2
:
rn1 = @reaction_network begin
hill(X,v,K,n), ∅ --> X
v*X^n/(X^n+K^n), ∅ --> X
end
\[ \begin{align*} \varnothing &\xrightarrow{\frac{v X^{n}}{K^{n} + X^{n}}} \mathrm{X} \\ \varnothing &\xrightarrow{\frac{v X^{n}}{K^{n} + X^{n}}} \mathrm{X} \end{align*} \]
rn2 = @reaction_network begin
mm(X,v,K), ∅ --> X
v*X/(X+K), ∅ --> X
end
\[ \begin{align*} \varnothing &\xrightarrow{\frac{v X}{K + X}} \mathrm{X} \\ \varnothing &\xrightarrow{\frac{X v}{K + X}} \mathrm{X} \end{align*} \]
Repressor Hill (hillr
) and Michaelis-Menten (mmr
) functions are also provided:
rn1 = @reaction_network begin
hillr(X,v,K,n), ∅ --> X
v*K^n/(X^n+K^n), ∅ --> X
end
\[ \begin{align*} \varnothing &\xrightarrow{\frac{v K^{n}}{K^{n} + X^{n}}} \mathrm{X} \\ \varnothing &\xrightarrow{\frac{v K^{n}}{K^{n} + X^{n}}} \mathrm{X} \end{align*} \]
rn2 = @reaction_network begin
mmr(X,v,K), ∅ --> X
v*K/(X+K), ∅ --> X
end
\[ \begin{align*} \varnothing &\xrightarrow{\frac{v K}{K + X}} \mathrm{X} \\ \varnothing &\xrightarrow{\frac{K v}{K + X}} \mathrm{X} \end{align*} \]
Please see the API Rate Laws section for more details.
Including non-species variables
Non-species state variables can be specified in the DSL using the @variables
macro. These are declared similarly to species. For example,
rn_with_volume = @reaction_network begin
@variables V(t)
k*V, 0 --> A
end
\[ \begin{align*} \varnothing &\xrightarrow{k V\left( t \right)} \mathrm{A} \end{align*} \]
creates a network with one species
species(rn_with_volume)
1-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
A(t)
and one non-species
nonspecies(rn_with_volume)
1-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
V(t)
giving two state variables, always internally ordered by species and then nonspecies:
states(rn_with_volume)
2-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
A(t)
V(t)
rn_with_volume
could then be extended with constraint equations for how V(t)
evolves in time, see the associated tutorial.
Specifying alternative time variables and/or extra independent variables
While the DSL defaults to allowing t
as the time variable, one can use the @ivs
macro to specify an alternative independent variable. For example, to make s
the default time variable one can say
rn_with_s = @reaction_network begin
@ivs s
@variables V(s)
@species B(s)
k, A + V*B --> C
end
Model ##ReactionSystem#4088
States (4):
B(s)
A(s)
C(s)
V(s)
Parameters (1):
k
where we see all states are now functions of s
.
Similarly, if one wants states to be functions of more than one independent variable, for example to encode a spatial problem, one can list more than one variable, i.e. @ivs t x y
. Here the first listed independent variable is always chosen to represent time. For example,
rn_with_many_ivs = @reaction_network begin
@ivs s x
@variables V1(s) V2(s,x)
@species A(s) B(s,x)
k, V1*A --> V2*B + C
end
Model ##ReactionSystem#4095
States (5):
A(s)
B(s, x)
C(s, x)
V1(s)
V2(s, x)
Parameters (1):
k
Here again s
will be the time variable, and any inferred species, C
in this case, are made functions of both variables, i.e. C(s, x)
.
Interpolation of Julia variables
The DSL allows Julia variables to be interpolated for the network name, within rate constant expressions, or for species/stoichiometry within reactions. Using the lower-level symbolic interface we can then define symbolic variables and parameters outside of the macro, which can then be used within expressions in the DSL (see the Programmatic Construction of Symbolic Reaction Systems tutorial for details on the lower-level symbolic interface). For example,
@parameters k α
@variables t
@species A(t)
spec = A
par = α
rate = k*A
name = :network
rn = @reaction_network $name begin
$rate*B, 2*$spec + $par*B --> $spec + C
end
\[ \begin{align*} 2 \mathrm{A} + \alpha \mathrm{B} &\xrightarrow{k A B} \mathrm{A} + \mathrm{C} \end{align*} \]
As the parameters k
and α
were pre-defined and appeared via interpolation, we did not need to declare them within the @reaction_network
macro, i.e. they are automatically detected as parameters:
parameters(rn)
2-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
k
α
as are the species coming from interpolated variables
species(rn)
3-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
B(t)
C(t)
A(t)
When using interpolation, expressions like 2$spec
won't work; the multiplication symbol must be explicitly included like 2*$spec
.