# FAQs

## How to index solution objects using symbolic variables and observables?

One can directly use symbolic variables to index into SciML solution objects. Moreover, observables can also be evaluated in this way. For example, consider the system

```
using Catalyst, DifferentialEquations, Plots
rn = @reaction_network ABtoC begin
(k₊,k₋), A + B <--> C
end
# initial condition and parameter values
setdefaults!(rn, [:A => 1.0, :B => 2.0, :C => 0.0, :k₊ => 1.0, :k₋ => 1.0])
```

Let's convert it to a system of ODEs, using the conservation laws of the system to eliminate two of the species:

`osys = convert(ODESystem, rn; remove_conserved = true)`

\[ \begin{align} \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} =& k_- \left( - A\left( t \right) + \Gamma_2 \right) - k_+ \left( A\left( t \right) + \Gamma_1 \right) A\left( t \right) \end{align} \]

Notice the resulting ODE system has just one ODE, while algebraic observables have been added for the two removed species (in terms of the conservation law constants, `Γ[1]`

and `Γ[2]`

)

`observed(osys)`

\[ \begin{align} B\left( t \right) =& A\left( t \right) + \Gamma_1 \\ C\left( t \right) =& - A\left( t \right) + \Gamma_2 \end{align} \]

Let's solve the system and see how to index the solution using our symbolic variables

```
oprob = ODEProblem(osys, [], (0.0, 10.0), [])
sol = solve(oprob, Tsit5())
```

Suppose we want to plot just species `C`

, without having to know its integer index in the state vector. We can do this using the symbolic variable `C`

, which we can get at in several ways

`sol[osys.C]`

or

```
@unpack C = osys
sol[C]
```

To evaluate `C`

at specific times and plot it we can just do

```
t = range(0.0, 10.0, length=101)
plot(t, sol(t, idxs = C), label = "C(t)", xlabel = "t")
```

If we want to get multiple variables we can just do

```
@unpack A, B = osys
sol(t, idxs = [A, B])
```

Plotting multiple variables using the SciML plot recipe can be achieved like

`plot(sol; idxs = [A, B])`

## How to disable rescaling of reaction rates in rate laws?

As explained in the Reaction rate laws used in simulations section, for a reaction such as `k, 2X --> 0`

, the generated rate law will rescale the rate constant, giving `k*X^2/2`

instead of `k*X^2`

for ODEs and `k*X*(X-1)/2`

instead of `k*X*(X-1)`

for jumps. This can be disabled when directly `convert`

ing a `ReactionSystem`

. If `rn`

is a generated `ReactionSystem`

, we can do

`osys = convert(ODESystem, rn; combinatoric_ratelaws=false)`

Disabling these rescalings should work for all conversions of `ReactionSystem`

s to other `ModelingToolkit.AbstractSystem`

s.

## How to use non-integer stoichiometric coefficients?

```
using Catalyst
rn = @reaction_network begin
k, 2.5*A --> 3*B
end
```

\[ \begin{align*} 2.5 \mathrm{A} &\xrightarrow{k} 3.0 \mathrm{B} \end{align*} \]

or directly via

```
@parameters k b
@variables t
@species A(t) B(t) C(t) D(t)
rx1 = Reaction(k,[B,C],[B,D], [2.5,1],[3.5, 2.5])
rx2 = Reaction(2*k, [B], [D], [1], [2.5])
rx3 = Reaction(2*k, [B], [D], [2.5], [2])
@named mixedsys = ReactionSystem([rx1, rx2, rx3], t, [A, B, C, D], [k, b])
osys = convert(ODESystem, mixedsys; combinatoric_ratelaws = false)
```

\[ \begin{align} \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} =& 0 \\ \frac{\mathrm{d} B\left( t \right)}{\mathrm{d}t} =& - 2 k B\left( t \right) - 5 \left( B\left( t \right) \right)^{2.5} k + \left( B\left( t \right) \right)^{2.5} k C\left( t \right) \\ \frac{\mathrm{d} C\left( t \right)}{\mathrm{d}t} =& - \left( B\left( t \right) \right)^{2.5} k C\left( t \right) \\ \frac{\mathrm{d} D\left( t \right)}{\mathrm{d}t} =& 5 k B\left( t \right) + 4 \left( B\left( t \right) \right)^{2.5} k + 2.5 \left( B\left( t \right) \right)^{2.5} k C\left( t \right) \end{align} \]

Note, when using `convert(ODESystem, mixedsys; combinatoric_ratelaws=false)`

the `combinatoric_ratelaws=false`

parameter must be passed. This is also true when calling `ODEProblem(mixedsys,...; combinatoric_ratelaws=false)`

. As described above, this disables Catalyst's standard rescaling of reaction rates when generating reaction rate laws, see also the Reaction rate laws used in simulations section. Leaving this keyword out for systems with floating point stoichiometry will give an error message.

For a more extensive documentation of using non-integer stoichiometric coefficients, please see the Symbolic Stochiometries section.

## How to set default values for initial conditions and parameters?

How to set defaults when using the `@reaction_network`

macro is described in more detail here. There are several ways to do this. Using the DSL, one can use the `@species`

and `@parameters`

options:

```
using Catalyst
sir = @reaction_network sir begin
@species S(t)=999.0 I(t)=1.0 R(t)=0.0
@parameters β=1e-4 ν=0.01
β, S + I --> 2I
ν, I --> R
end
```

```
Model sir
States (3):
S(t) [defaults to 999.0]
I(t) [defaults to 1.0]
R(t) [defaults to 0.0]
Parameters (2):
β [defaults to 0.0001]
ν [defaults to 0.01]
```

When directly constructing a `ReactionSystem`

, we can set the symbolic values to have the desired default values, and this will automatically be propagated through to the equation solvers:

```
using Catalyst, Plots, OrdinaryDiffEq
@parameters β=1e-4 ν=.01
@variables t
@species S(t)=999.0 I(t)=1.0 R(t)=0.0
rx1 = Reaction(β, [S, I], [I], [1,1], [2])
rx2 = Reaction(ν, [I], [R])
@named sir = ReactionSystem([rx1, rx2], t)
oprob = ODEProblem(sir, [], (0.0, 250.0))
sol = solve(oprob, Tsit5())
plot(sol)
```

One can also build a mapping from symbolic parameter/species to value/initial condition and pass these to the `ReactionSystem`

via the `defaults`

keyword argument:

```
@parameters β ν
@variables t
@species S(t) I(t) R(t)
rx1 = Reaction(β, [S,I], [I], [1,1], [2])
rx2 = Reaction(ν, [I], [R])
defs = [β => 1e-4, ν => .01, S => 999.0, I => 1.0, R => 0.0]
@named sir = ReactionSystem([rx1, rx2], t; defaults = defs)
```

Finally, default values can also be added after creating the system via the `setdefaults!`

command and passing a `Symbol`

based mapping, like

```
sir = @reaction_network sir begin
β, S + I --> 2I
ν, I --> R
end
setdefaults!(sir, [:β => 1e-4, :ν => .01, :S => 999.0, :I => 1.0, :R => 0.0])
```

## How to specify initial conditions and parameters values for `ODEProblem`

and other problem types?

To explicitly pass initial conditions and parameters we can use mappings from Julia `Symbol`

s corresponding to each variable/parameter to their values, or from ModelingToolkit symbolic variables/parameters to their values. Using `Symbol`

s we have

```
using Catalyst, DifferentialEquations
rn = @reaction_network begin
α, S + I --> 2I
β, I --> R
end
u0 = [:S => 999.0, :I => 1.0, :R => 0.0]
p = (:α => 1e-4, :β => .01)
op1 = ODEProblem(rn, u0, (0.0, 250.0), p)
```

while using ModelingToolkit symbolic variables we have

```
@parameters α β
@variables t
@species S(t) I(t) R(t)
u0 = [S => 999.0, I => 1.0, R => 0.0]
p = (α => 1e-4, β => .01)
op2 = ODEProblem(rn, u0, (0.0, 250.0), p)
```

*Note,* while symbolic mappings as in the last example will work with *any* `ModelingToolkit.AbstractSystem`

, for example if one `convert`

s `rn`

to an `ODESystem`

, `Symbol`

-based mappings only work when passing a `ReactionSystem`

directly into a problem type. That is, the following does not work

```
osys = convert(ODESystem, rn)
# this fails
u0 = [:S => 999.0, :I => 1.0, :R => 0.0]
p = (:α => 1e-4, :β => .01)
op = ODEProblem(osys, u0, (0.0, 250.0), p)
```

In this case one must either use a symbolic mapping as was used to make `op2`

in the second example, or one can use the `symmap_to_varmap`

function to convert the `Symbol`

mapping to a symbolic mapping. I.e. this works

```
osys = convert(ODESystem, rn)
# this works
u0 = symmap_to_varmap(rn, [:S => 999.0, :I => 1.0, :R => 0.0])
p = symmap_to_varmap(rn, (:α => 1e-4, :β => .01))
op = ODEProblem(osys, u0, (0.0, 250.0), p)
```

## How to include non-reaction terms in equations for a chemical species?

One method to add non-reaction terms into an ODE or algebraic equation for a chemical species is to add a new (non-species) state variable that represents those terms, let it be the rate of zero order reaction, and add a constraint equation. I.e., to add a force of `(1 + sin(t))`

to $dA/dt$ in a system with the reaction `k, A --> 0`

, we can do

```
using Catalyst
@variables t f(t)
rx1 = @reaction k, A --> 0
rx2 = @reaction $f, 0 --> A
eq = f ~ (1 + sin(t))
@named rs = ReactionSystem([rx1, rx2, eq], t)
osys = convert(ODESystem, rs)
```

\[ \begin{align} \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} =& f\left( t \right) - k A\left( t \right) \\ f\left( t \right) =& 1 + \sin\left( t \right) \end{align} \]

In the final ODE model, `f`

can be eliminated by using `ModelingToolkit.structural_simplify`

```
osyss = structural_simplify(osys)
full_equations(osyss)
```

\[ \begin{align} \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} =& 1 + \sin\left( t \right) - k A\left( t \right) \end{align} \]

## How to modify generated ODEs?

Conversion to other `ModelingToolkit.AbstractSystem`

s allows the possibility to modify the system with further terms that are difficult to encode as a chemical reaction or a constraint equation. For example, an alternative method to the previous question for adding a forcing function, $1 + \sin(t)$, to the ODE for `dA/dt`

is

```
using Catalyst
rn = @reaction_network begin
k, A --> 0
end
osys = convert(ODESystem, rn)
dAdteq = equations(osys)[1]
t = ModelingToolkit.get_iv(osys)
dAdteq = Equation(dAdteq.lhs, dAdteq.rhs + 1 + sin(t))
# create a new ODESystem with the modified equation
@named osys2 = ODESystem([dAdteq], t)
```

\[ \begin{align} \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} =& 1 + \sin\left( t \right) - k A\left( t \right) \end{align} \]

## How to override mass action kinetics rate laws?

While generally one wants the reaction rate law to use the law of mass action, so the reaction

```
using Catalyst
rn = @reaction_network begin
k, X --> ∅
end
convert(ODESystem, rn)
```

\[ \begin{align} \frac{\mathrm{d} X\left( t \right)}{\mathrm{d}t} =& - k X\left( t \right) \end{align} \]

occurs at the (ODE) rate $d[X]/dt = -k[X]$, it is possible to override this by using any of the following non-filled arrows when declaring the reaction: `<=`

, `⇐`

, `⟽`

, `=>`

, `⇒`

, `⟾`

, `⇔`

, `⟺`

. This means that the reaction

```
rn = @reaction_network begin
k, X => ∅
end
convert(ODESystem, rn)
```

\[ \begin{align} \frac{\mathrm{d} X\left( t \right)}{\mathrm{d}t} =& - k \end{align} \]

will occur at rate $d[X]/dt = -k$ (which might become a problem since $[X]$ will be degraded at a constant rate even when very small or equal to 0).

Note, stoichiometric coefficients are still included, i.e. the reaction

```
rn = @reaction_network begin
k, 2*X ⇒ ∅
end
convert(ODESystem, rn)
```

\[ \begin{align} \frac{\mathrm{d} X\left( t \right)}{\mathrm{d}t} =& - 2 k \end{align} \]

has rate $d[X]/dt = -2 k$.

## How to specify user-defined functions as reaction rates?

The reaction network DSL can "see" user-defined functions that work with ModelingToolkit. e.g., this is should work

```
using Catalyst
myHill(x) = 2*x^3/(x^3+1.5^3)
rn = @reaction_network begin
myHill(X), ∅ → X
end
```

\[ \begin{align*} \varnothing &\xrightarrow{\frac{2 X^{3}}{3.375 + X^{3}}} \mathrm{X} \end{align*} \]

In some cases, it may be necessary or desirable to register functions with Symbolics.jl before their use in Catalyst, see the discussion here.