Parametric Stoichiometry

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

Note, this tutorial requires ModelingToolkit v8.5.4 or greater to work properly.

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 k₊ k₋ m n
2-element Vector{Reaction}:
 k₊, m*A --> (m*n)*B
 k₋, B --> A

Note, as always the @reaction_network macro sets all symbols not declared to be parameters to be species, so that in this example we have two species, A and B, and four parameters. In addition, the stoichiometry is applied to the right most 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 k
1-element Vector{Reaction}:
 k, (A(t))*C --> 2*B

we see three species, (A,B,C), however, A is treated as the stoichiometric coefficient of C, i.e.

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

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, A(t), B(t)
rxs = [Reaction(k₊,[A],[B],[m],[m*n]),
revsys2 = ReactionSystem(rxs,t; name=:revsys)
revsys2 == revsys

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

Note, the @reaction macro assumes all symbols are parameters except the right most symbols in the reaction line (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)
2-element Vector{Equation}:
 Differential(t)(A(t)) ~ k₋*B(t) + (-k₊*m*(A(t)^m)) / factorial(m)
 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)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1.0)
u0: 2-element Vector{Float64}:

We can now solve and plot the system

sol = solve(oprob, Tsit5())

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)
2-element Vector{Equation}:
 Differential(t)(A(t)) ~ k₋*B(t) - k₊*m*(A(t)^m)
 Differential(t)(B(t)) ~ k₊*m*n*(A(t)^m) - k₋*B(t)

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())

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 k₊ k₋ kₚ γₚ
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)
(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.Mul{Real, Int64, Dict{Any, Number}, Nothing}, Vector{Equation}}(k₋*(P(t)^2)*G₊(t), Equation[G₊(t) ~ G₊(t) - 1, G₋(t) ~ 1 + G₋(t)]), JumpProcesses.ConstantRateJump{SymbolicUtils.Mul{Real, Int64, Dict{Any, Number}, Nothing}, Vector{Equation}}(kₚ*G₊(t), Equation[P(t) ~ 1 + P(t) + rand(Distributions.Geometric(1 / b))])], 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

1-element Vector{Equation}:
 P(t) ~ 1 + P(t) + rand(Distributions.Geometric(1 / b))

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)
    Pmean ./= Nsims
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).