# Symbolic Stochiometries

Catalyst supports stoichiometric coefficients that involve parameters, species, or even general expressions. In this tutorial we show several examples of how to use symbolic stoichiometries, and discuss several caveats to be aware of.

## Using symbolic stoichiometry

Let's first consider a simple reversible reaction where the number of reactants is a parameter, and the number of products is the product of two parameters.

using Catalyst, Latexify, DifferentialEquations, ModelingToolkit, Plots
revsys = @reaction_network revsys begin
k₊, m*A --> (m*n)*B
k₋, B --> A
end
reactions(revsys)
2-element Vector{Reaction}:
k₊, m*A --> (m*n)*B
k₋, B --> A

Note, as always the @reaction_network macro defaults to setting all symbols neither used as a reaction substrate nor a product to be parameters. Hence, in this example we have two species (A and B) and four parameters (k₊, k₋, m, and n). In addition, the stoichiometry is applied to the rightmost symbol in a given term, i.e. in the first equation the substrate A has stoichiometry m and the product B has stoichiometry m*n. For example, in

rn = @reaction_network begin
k, A*C --> 2B
end
reactions(rn)
1-element Vector{Reaction}:
k, A*C --> 2*B

we see two species, (B,C), with A treated as a parameter representing the stoichiometric coefficient of C, i.e.

rx = reactions(rn)[1]
rx.substrates[1],rx.substoich[1]
(C(t), A)

We could have equivalently specified our systems directly via the Catalyst API. For example, for revsys we would could use

@parameters k₊, k₋, m, n
@variables t
@species A(t), B(t)
rxs = [Reaction(k₊, [A], [B], [m], [m*n]),
Reaction(k₋, [B], [A])]
revsys2 = ReactionSystem(rxs,t; name=:revsys)
revsys2 == revsys
true

which can be simplified using the @reaction macro to

rxs2 = [(@reaction k₊, m*A --> (m*n)*B),
(@reaction k₋, B --> A)]
revsys3 = ReactionSystem(rxs2,t; name=:revsys)
revsys3 == revsys
true

Note, the @reaction macro again assumes all symbols are parameters except the substrates or reactants (i.e. A and B). For example, in @reaction k, F*A + 2(H*G+B) --> D, the substrates are (A,G,B) with stoichiometries (F,2*H,2).

Let's now convert revsys to ODEs and look at the resulting equations:

osys = convert(ODESystem, revsys)
equations(osys)
2-element Vector{Equation}:
Differential(t)(A(t)) ~ (-k₊*m*(A(t)^m)) / factorial(m) + k₋*B(t)
Differential(t)(B(t)) ~ (k₊*m*n*(A(t)^m)) / factorial(m) - k₋*B(t)

Notice, as described in the Reaction rate laws used in simulations section, the default rate laws involve factorials in the stoichiometric coefficients. For this reason we must specify m and n as integers, and hence use a tuple for the parameter mapping

p  = (k₊ => 1.0, k₋ => 1.0, m => 2, n => 2)
u₀ = [A => 1.0, B => 1.0]
oprob = ODEProblem(osys, u₀, (0.0, 1.0), p)

We can now solve and plot the system

sol = solve(oprob, Tsit5())
plot(sol)

If we had used a vector to store parameters, m and n would be converted to floating point giving an error when solving the system.

An alternative approach to avoid the issues of using mixed floating point and integer variables is to disable the rescaling of rate laws as described in Reaction rate laws used in simulations section. This requires passing the combinatoric_ratelaws=false keyword to convert or to ODEProblem (if directly building the problem from a ReactionSystem instead of first converting to an ODESystem). For the previous example this gives the following (different) system of ODEs

osys = convert(ODESystem, revsys; combinatoric_ratelaws = false)
equations(osys)
2-element Vector{Equation}:
Differential(t)(A(t)) ~ k₋*B(t) - k₊*m*(A(t)^m)
Differential(t)(B(t)) ~ -k₋*B(t) + k₊*m*n*(A(t)^m)

Since we no longer have factorial functions appearing, our example will now run even with floating point values for m and n:

p  = (k₊ => 1.0, k₋ => 1.0, m => 2.0, n => 2.0)
oprob = ODEProblem(osys, u₀, (0.0, 1.0), p)
sol = solve(oprob, Tsit5())
plot(sol)

## Gene expression with randomly produced amounts of protein

As a second example, let's build the negative feedback model from MomentClosure.jl that involves a bursty reaction that produces a random amount of protein.

In our model G₋ will denote the repressed state, and G₊ the active state where the gene can transcribe. P will denote the protein product of the gene. We will assume that proteins are produced in bursts that produce m proteins, where m is a (shifted) geometric random variable with mean b. To define m we must register the Distributions.Geometric distribution from Distributions.jl with Symbolics.jl, after which we can use it in symbolic expressions:

using Distributions: Geometric
@register_symbolic Geometric(b)
@parameters b
m = rand(Geometric(1/b)) + 1

Note, as we require the shifted geometric distribution, we add one to Distributions.jl's Geometric random variable (which includes zero).

We can now define our model

burstyrn = @reaction_network burstyrn begin
k₊, G₋ --> G₊
k₋*P^2, G₊ --> G₋
kₚ, G₊ --> G₊ + \$m*P
γₚ, P --> ∅
end
reactions(burstyrn)
4-element Vector{Reaction}:
k₊, G₋ --> G₊
k₋*(P(t)^2), G₊ --> G₋
kₚ, G₊ --> G₊ + (1 + rand(Distributions.Geometric(1 / b)))*P
γₚ, P --> ∅

The parameter b does not need to be explicitly declared in the @reaction_network macro as it is detected when the expression rand(Geometric(1/b)) + 1 is substituted for m.

We next convert our network to a jump process representation

jsys = convert(JumpSystem, burstyrn; combinatoric_ratelaws = false)
equations(jsys)
(JumpProcesses.MassActionJump[JumpProcesses.MassActionJump{Num, Vector{Pair{Any, Int64}}, Vector{Pair{Any, Int64}}, Nothing}(k₊, Pair{Any, Int64}[G₋(t) => 1], Pair{Any, Int64}[G₊(t) => 1, G₋(t) => -1], nothing), JumpProcesses.MassActionJump{Num, Vector{Pair{Any, Int64}}, Vector{Pair{Any, Int64}}, Nothing}(γₚ, Pair{Any, Int64}[P(t) => 1], Pair{Any, Int64}[P(t) => -1], nothing)], JumpProcesses.ConstantRateJump[JumpProcesses.ConstantRateJump{SymbolicUtils.BasicSymbolic{Real}, Vector{Equation}}(k₋*G₊(t)*(P(t)^2), Equation[G₊(t) ~ -1 + G₊(t), G₋(t) ~ 1 + G₋(t)]), JumpProcesses.ConstantRateJump{SymbolicUtils.BasicSymbolic{Real}, Vector{Equation}}(kₚ*G₊(t), Equation[P(t) ~ 1 + rand(Distributions.Geometric(1 / b)) + P(t)])], JumpProcesses.VariableRateJump[])

Notice, the equations of jsys have three MassActionJumps for the first three reactions, and one ConstantRateJump for the last reaction. If we examine the ConstantRateJump more closely we can see the generated rate and affect! functions for the bursty reaction that makes protein

equations(jsys)[4].rate
kₚ*G₊(t)
equations(jsys)[4].affect!
1-element Vector{Equation}:
P(t) ~ 1 + rand(Distributions.Geometric(1 / b)) + P(t)

Finally, we can now simulate our jumpsystem

pmean = 200
bval = 70
γₚval = 1
k₋val = 0.001
k₊val = 0.05
kₚval = pmean * γₚval * (k₋val * pmean^2 + k₊val) / (k₊val * bval)
p = symmap_to_varmap(jsys, (:k₊ => k₊val, :k₋ => k₋val, :kₚ => kₚval,
:γₚ => γₚval, :b => bval))
u₀ = symmap_to_varmap(jsys, [:G₊ => 1, :G₋ => 0, :P => 1])
tspan = (0., 6.0)   # time interval to solve over
dprob = DiscreteProblem(jsys, u₀, tspan, p)
jprob = JumpProblem(jsys, dprob, Direct())
sol = solve(jprob, SSAStepper())
plot(sol.t, sol[jsys.P], legend = false, xlabel = "time", ylabel = "P(t)")

To double check our results are consistent with MomentClosure.jl, let's calculate and plot the average amount of protein (which is also plotted in the MomentClosure.jl tutorial).

function getmean(jprob, Nsims, tv)
Pmean = zeros(length(tv))
@variables t, P(t)
for n in 1:Nsims
sol = solve(jprob, SSAStepper())
Pmean .+= sol(tv, idxs=P)
end
Pmean ./= Nsims
end
tv = range(tspan[1],tspan[2],step=.1)
psim_mean = getmean(jprob, 20000, tv)
plot(tv, psim_mean; ylabel = "average of P(t)", xlabel = "time",
xlim = (0.0,6.0), legend = false)

Comparing, we see similar averages for P(t).